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ABSTRACT 


This  thesis  describes  the  development  of  a  computer  program  for 
analysing  the  flow  of  water  through  soils  which  consist  of  a  number  of  zones 
having  different  permeability;  the  method  of  analysis  provides  for  each  zone 
being  anisotropic  to  a  different  degree,  and  for  the  directions  of  principal 
permeability  being  different  in  each  zone. 

The  analysis  involves  finding  a  numerically  approximate  solution  of 
the  basic  equation  of  flow  at  every  node  in  a  mesh  superimposed  on  the  flow 
region.  To  avoid  difficulties  which  occur  in  the  conventional  iteration  pro¬ 
cedure  at  zone  boundaries,  the  mesh  has  three  axes;  each  axis  is  parallel  to 
the  sides  of  an  irregular  triangle,  drawn  so  that  it  coincides  with  zone 
boundaries . 


The  finite  difference  form  of  the  basic  flow  equation  in  terms  of  three 
axes  in  a  plane  is  derived  from  first  principles.  The  assembly  of  the  computer 
program  is  fully  described,  and  the  results  of  several  problems  for  which  the 
program  has  been  used  are  presented.  As  presently  developed,  the  program 
applies  to  steady-state,  two-dimensional,  confined  flow,  but  it  is  suggested 
that  it  could  be  methodically  extended  to  cover  cases  outside  these  restrictions. 

The  writer  submits  that  the  method  enables  analyses  to  be  made  of 
flow  regions  which  are  usually  regarded  as  too  complicated  for  solution  by 
existing  methods,  unless  simplifying  approximations  are  made. 
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GL05SARY  OF  SYMBOLS 


A  list  of  the  symbols  used  in  the  computer  program  is  given 
on  pages  D2  to  D4. 
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Coefficients  to  be  used  in  the  finite  difference  equations. 

A  second  subscript  signifies  the  number  of  the  triangle  with 
which  the  coefficient  is  associated. 


f  Superscript  indicating  that  a  value  is  ficticious,  i.e.  occurs 

at  an  image  point. 

P  The  ratio  of  a  side  of  a  triangle  to  the  sine  of  an  opposite 

angle. 

h  Total  head  of  water. 


h 
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Total  head  of  water  at  a  node  under  consideration. 


h' 
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h' 
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h" 
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h" 
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h" 
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Potential  heads  at  adjacent  nodes  to  the  node  under  consideration 
in  positive  and  negative  directions  along  the  Ou,  Ov  and  Ow  axes. 
A  second  subscript  signifies  the  number  of  the  triangle  in  which 
the  node  is  situated. 


(h) 


u 


The  partial  derivative  of  h  with  respect  to  u. 


(h)  The  second  partial  derivative  of  h  with  respect  to  u  and  v. 


k  Permeability. 

M  p  Symbols  substituted  for  u,  v,  w,  and  vf  w,  u,  and 

*  *  w,  u,  v,  respectively. 

q  Flow  through  an  elemental  arc  of  a  circle  radius  As. 

Outward  flow  across  the  curved  portion  of  the  perimeter  of  a 
1  semi-circle  radius  As* 


( viii ) 


. 


(ix) 
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U,  V,  w 


V 


x»  y 


xi  t  yi  » 

zi 


Inward  flow  across  the  diameter  of  a  semi-circle  radius  As. 


Residual  at  a  node. 


Any  general  direction  from  the  centre  of  an  elemental  circle. 


Symbols  signifying  each  of  the  three  axes  in  a  triangular  mesh. 


Superficial  velocity. 


Directions  of  principal  permeability  in  an  anisotropic  soil. 


Three  mutually  perpendicular  directions  in  an  isotropic  soil. 


As 


Radius  of  an  elemental  circle. 


A-|U 
A-j  v 
A)W 
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•Node  spacings  on  the  three  axes  of  a  triangular  mesh. 


Angle  between  the  direction  of  minimum  permeability  and  any  general 
direction,  s. 


<*a 


Angle  between  the  horizontal  and  the  direction  of  minimum  permeability. 


°<u 

<Xv 


«w 


Angles  between  the  direction  of  minimum  permeability 
Ov  and  Ow  axes. 


and  the  Ou, 


,  Angles  between  the  horizontal  and  the  Ou,  Ov  and  Ow  axes. 


CHAPTER  I 


INTRODUCTION 


1.1  General 

The  flow  of  water  through  a  soil  may  be  analysed  by  one  of  several 
methods.  The  method  chosen  is  governed  largely  by  the  accuracy  required. 

The  hydraulic  conditions  set  up  when  water  flows  through  soil  under  ideal 
circumstances  can  be  described  concisely  by  simple  equations;  consequently 
flow  analyses  are  amenable  to  mathematical  treatment.  However,  estimates  of 
permeabilities  of  different  soils  at  a  site  are  frequently  only  crude  and  thus 
an  exact  mathematical  analysis  is  seldom  warranted;  also,  the  boundary  con¬ 
ditions  of  most  practical  problems  are  sufficiently  complex  to  make  the  analy¬ 
sis  very  difficult  by  classical  methods  of  mathematics. 

The  engineer  is  then  forced  into  making  decisions  such  as:  he  will 
accept  a  less  rigorous  method  of  analysis  in  the  light  of  the  uncertainty  in 
the  estimates  of  permeability  and  the  boundary  conditions  of  the  soil;  he  will 
treat  a  certain  soil  as  infinitely  permeable  compared  with  another  soil  at  the 
same  site  if  the  ratios  of  permeability  of  the  two  are  more  than,  say,  ten 
to  one;  he  will  treat  a  soil  as  if  it  were  isotropic  even  though  it  may  be 
anisotropic.  In  many  instances,  the  effect  of  these  decisions  on  the  final 
result  will  be  sufficiently  small  to  be  negligible,  especially  if  the  quan¬ 
tity  of  seepage  is  all  that  is  required.  If  estimates  of  critical  hydraulic 
gradients  or  distribution  of  pore  pressures  are  required,  more  careful  thought 
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may  be  taken  in  choosing  a  method  of  analysis  and  in  making  any  simplifying 
assumptions • 

To  reach  a  practical  solution,  the  engineer  almost  always  has  to 
estimate  values  for  some  missing  soil  data,  and  modify  the  data  that  is 
available,  so  that  the  final  answer  is  a  mixture  of  intuition  and  analysis. 

The  writer  has  endeavoured  to  devise  a  method  of  solution  which  reduces  the 
intuition  factor,  and  yet  is  still  practical, 

1,2  Statement  of  the  problem 

In  view  of  the  inroads  that  the  digital  computer  is  making  on  engin¬ 
eering  science  at  the  present  time,  the  writer  became  inclined  to  investigate 
its  use  in  regard  to  seepage  problems, 

A  method  of  analysis  is  available  which  is  adaptable  to  the  computer 
(Scott,  1963,  page  133);  this  method  involves  a  repetitive  process  of  calcu¬ 
lations  and  the  end  result  is  approached  through  successively  better  approxi¬ 
mations,  If  the  boundary  conditions  of  the  flow  region  are  simple  (i.e,,  in 
the  extreme,  if  the  flow  region  is  rectangular)  and  if  the  soil  is  homogeneous 
and  isotropic,  the  length  of  the  computer  program  using  this  method  is  small. 
However,  for  such  a  case,. the  result  could  be  achieved  quite  easily  by  hand, 
using  the  same  method;  alternatively  another  method  of  analysis  could  be  used 
to  obtain  an  equally  good  result.  If  the  flow  region  consists  of  a  number  of 
soils  of  different  permeabilities,  as  may  be  found  in  a  zoned  earth  dam  or  a 
natural  soil  deposit  with  several  layers  of  different  materials,  the  method 
is  much  more  complex  and  the  assembling  of  the  anpropriate  computer  program 
is  time-consuming.  Moreover,  if  another  flow  region  of  similar  complexity  is 
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then  considered,  the  time  spent  in  altering  the  program  to  fit  the  con¬ 
ditions  of  the  second  case  may  make  it  worthwhile  to  assemble  an  entirely 
new  program.  The  approach  to  the  problem  using  other  methods  would  be 
extremely  difficult  or  impossible,  unless  simplifying  assumptions  were 
made . 

The  problem  that  confronted  the  writer,  and  to  which  this  thesis 
is  devoted,  was  to  assemble  a  computer  program  which  would  be  applicable 
to  seepage  through  a  non-homogeneous , ^  anisotropic  soil  section  of  any 
desired  configuration  and  complexity.  The  program  could  then  be  used  for  all 
types  of  flow  regions,  without  any  alterations. 

1.3  Extent  and  limitations  of  work  carried  out 

The  problem  was  approached  by  modifying  the  conventional  procedures 
involved  in  the  repetitive  process  of  calculations  so  that  cumbersome  steps 
caused  by  the  presence  of  irregular  boundaries  would  be  eliminated.  In  this 
way,  the  number  and  configuration  of  different  zones  of  soil  was  a  factor  of 
minor  importance  in  the  programming,  instead  of,  as  hitherto,  a  major  difficulty 
The  conventional  procedure  uses  the  basic  flow  equation  in  two  dimensions;  this 
equation  is  used  herein  but  in  a  rather  unusual  form.  One  outcome  of  this  modif 
ication  is  that,  if  the  solution  is  attempted  by  hand,  even  for  the  simple  case 
of  the  rectangular  flow  region  mentioned  in  SECTION  1.2,  the  calculations  be¬ 
come  very  laborious.  Laborious  calculations  present  no  problems  in  programming 
and  computer  work,  but  cumbersome  operations  are  to  be  avoided  if  possible, 

soil  which  consists  of  two  or  more  zones,  homogeneous  in  themselves,  but 
having  different  permeabilities  from  each  other,  is  regarded  as  non-homogeneous 
and  this  meaning  is  used  throughout  this  thesis . 
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The  computer  program  was  assembled  so  that  the  total  heads  at 
reasonably  closely  spaced  points  within  the  flow  region  would  appear  as 
output,  together  with  coordinates  to  locate  the  positions  of  the  points. 
The  program  applies  to  two-dimensional  steady-state  flow  through  a  non- 
homogeneous,  anisotropic  soil  for  which: 

(a)  Darcy's  law  and  continuity  conditions  apply  and  thus  the  flow 
may  be  described  by  the  basic  equation: 
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where  x  and  y  denote  the  directions  of  principal  permeabilities  in 
an  anisotropic  soil, 

k*>y  are  the  coefficients  of  permeability  in  the  x-  and  y- 
directions  respectively, 

h  is  the  total  head  of  water  at  an  arbitary  point  within 
the  mass, 

(b)  inside  the  boundaries  of  any  zone,  the  soil  is  homogeneous,  but  may 
be  anisotropic, 


(c)  the  boundaries  of  each  zone  are  defined, 

(d)  for  an  anisotropic  soil,  the  principal  coefficients  of  permeability 
and  the  direction  in  which  they  act  are  defined;  for  an  isotropic  soil,  the 
coefficient  of  permeability  is  defined, 

(e)  none  of  the  flow  region  boundaries  are  controlled  by  gravity,  i.e., 
there  is  no  phreatic  line  present, 

(f)  the  flow  region  boundaries  are  finite  in  length  and  are  defined. 


Item  (a)  above  always  applies  in  practical  seepage  problems.  Darcy's 

law  is  known  to  hold  if  the  seepage  velocity  of  the  water  is  such  that  the 

2 

Reynold's  number  is  less  than  about  one;  continuity  conditions  are  upheld  if 
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All  symbols  have  been  listed  in  the  Glossary  of  Symbols. 

The  Reynold's  number  in  this  context  is  based  on  superficial  velocity  and  a  diameter 
which  i3  that  of  a  sphere  having  a  volume  equal  to  the  quotient  of  the  volume  of 
solids  of  the  sample  and  the  number  of  grains  in  the  sample. 
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the  soil  is  saturated  within  the  flow  region,  since  water  and  soil  can  be 
considered  incompressible  for  practical  cases.  In  order  to  satisfy  item 
(b)  for  cases  in  which  the  permeability  varies  gradually  from  point  to  point 
within  a  zone,  the  zone  should  be  sub-divided  into  a  number  of  smaller  zones. 
Items  (c)  and  (d)  are  frequently  incomplete  as  a  result  of  a  site  investi¬ 
gation  alone.  The  engineer  will  usually  have  to  make  up  the  deficiencies 
with  estimated  data  obtained  by  informed  guesswork;  he  may  find  it  profitable 
to  run  the  program  through  the  computer  several  times  for  the  same  problem, 
each  time  inserting  different  values  for  the  guessed  data;  by  comparing  the 
predicted  flow  patterns  for  each  set  of  data,  the  worst  conditions  may  be 
determined.  Item  (e)  is  a  limitation  caused  by  lack  of  time  available  to  the 
writer.  The  main  effect  of  this  limitation  is  that  the  program,  as  it  stands, 
cannot  be  used  for  analysing  the  flow  of  water  through  earth  dams.  The  writer 
believes  that  the  program  could  be  extended,  without  much  difficulty,  so  that 
flow  problems  containing  a  phreatic  line  may  be  solved;  such  an  extension  need 
not  involve  any  loss  of  generality.  The  effect  of  item  (f)  is  that  the  flow 
region  must  be  enclosed  within  a  finite  area.  In  a  strict  sense,  such  is  not 
the  case  when  the  problem  involves  flow  under  a  spillway  situated  on  an  infin¬ 
ite  layer  of  permeable  material  overlying  a  horizontal  impermeable  base.  To 
analyse  such  a  problem,  vertical  impermeable  boundaries  would  have  to  be  in¬ 
serted  through  the  permeable  layer  at  some  finite  distance  from  the  spillway; 
if  the  distance  were  made  sufficiently  large,  the  effect  of  the  presence  of 
these  boundaries  would  be  negligible. 

Limitations  to  the  size  of  the  data  for  any  flow  region  are  listed 
below.  These  limitations  were  fixed  by  the  number  of  memory  storage  locations 
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available  in  the  computer  used  by  the  writer.  For  most  practical  problems, 
these  limitations  will  have  no  effect,  but  should  a  problem  be  encountered 
where  the  size  of  the  data  is  greater  than  that  specified,  the  program  may 
be  modified  so  that  the  computer  calls  upon  magnetic  tape  for  storing  the 
additional  information.  The  limitations  are  that; 

(a)  the  maximum  number  of  zones  is  20, 

(b)  considering  the  shape  of  any  zone  to  be  a  polygon,  the  maximum 
number  of  boundaries  of  any  zone  is  20, 

(c)  the  maximum  number  of  impermeable  boundaries  of  the  flow  region  is 
30, 

(d)  the  maximum  number  of  boundaries  across  which  water  enters  or  leaves 
the  flow  region  (i.e.-,  constant  total  head  boundaries)  is  30. 

In  the  use  of  the  program,  any  curved  boundaries  which  exist  in  the 
flow  region  must  be  approximated  by  a  series  of  straight  lines.  Since  the  data 
defining  zone  boundaries  is  presented  to  the  computer  by  simply  listing  coordin¬ 
ates  locating  the  corner  points  of  the  polygons,  it  does  not  matter  how  irreg¬ 
ular  the  polygons  are.  The  ratios  of  principal  permeabilities  for  the  soils 
in  each  zone  may  be  different,  and  the  direction  of  stratification  in  an  aniso¬ 
tropic  soil  in  any  zone  may  be  different  from  that  in  any  other  zone.  These 
conditions  may  hardly  ever  be  realized  in  practice,  but  the  extra  labour  involved 
in  programming  for  them  is  minimal  and  if  they  were  not  included  in  the  program, 
some  generality  would  be  lost. 

A  review  of  present  methods  has  been  undertaken;  the  method  of  analy¬ 
sis  developed  by  the  writer  is  described  herein,  together  with  the  computer  pro¬ 
gram.  Some  examples  for  which  the  program  has  been  used  are  presented. 


CHAPTER  II 


REVIEW  OF  EXISTING  METHODS  OF  FLOW  ANALYSIS 


2.1  Basis  of  flow  analyses 

Methods  for  analysing  seepage  problems  are  based  on  the  condition 
that,  when  water  flows  through  a  homogeneous,  isotropic  soil  mass,  the  Laplace 
equation,  set  out  as  follows,  is  satisfied: 


atf  +  a*,* 


-  O 


(2.1a) 


where  h  is  the  total  head  of  water  an  an  arbitary  point  within  the  mass, 
and  x^y^z^  are  three  mutually  perpendicular  directions. 

The  derivation  of  equation  (2.1a)  depends  on  the  assumption  that  Darcy's  law 
and  continuity  conditions  are  upheld.  In  many  practical  cases  in  which  the 
length  of  the  flow  region  in  a  direction  transverse  to  the  flow  is  long  com¬ 
pared  with  the  width  or  depth,  the  flow  conditions  may  be  considered  as  two- 
dimensional.  Equation  (2.1a)  for  such  cases  becomes  modified  to: 


<£h  ^  £h 


O 


(2.1b) 


where 


and  y^  are  perpendicular  directions  in  the  plane  parallel 
to  flow. 


In  an  isotropic  soil  through  which  water  is  percolating,  flow  occurs 
at  any  point  in  the  direction  of  the  maximum  hydraulic  gradient  at  that  point. 
After  obtaining  the  solution  to  equation  (2.1b),  curves  may  be  plotted  joining 
points  at  which  the  total  heads  are  the  same;  the  maximum  hydraulic  gradient 
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at  any  point  acts  in  a  direction  normal  to  the  curve  of  constant  head 
passing  through  that  point*  Hence  curves  which  represent  the  path  of  an 
individual  water  particle  may  be  superimposed  on  several  curves  of  constant 
head  by  ensuring  that  they  intersect  each  other  at  right  angles.  The  curves 
of  constant  head  and  those  representing  the  flow  paths  are  called,  respectively, 
equipotential  lines  and  flow  lines;  the  combined  representation  of  equipotential 
lines  and  flow  lines  is  called  a  flow  net. 

2.2  Methods  for  analysing  flow 

The  available  methods  for  solving  flow  problems  may  be  broadly  divi¬ 
ded  into  three  catagories: 

(a)  Theoretical  methods 

These  methods  entail  solving,  either  directly  or  by  numerical  approxi¬ 
mation,  differential  equations  representing  the  equipotential  lines  and  flow 
lines.  The  direct  method  uses  the  theory  of  complex  variables  to  study  the 
orthogonality  of  the  two  sets  of  curves.  Briefly,  by  choosing  an  appropriate 
equation,  the  geometry  of  the  flow  region  can  be  transformed  into  another  con¬ 
figuration,  or  map,  the  flow  net  in  which  is  known  or  is  obvious.  Each  point 
in  the  flow  region  before  transformation  is  represented  by  one  point  in  the 
mapped  region,  and,  therefore,  having  ascertained  the  flow  conditions  at  any 
point  in  the  mapped  region,  these  conditions  can  be  related  to  the  corresponding 
point  in  the  original  flow  region.  The  application  of  the  complex  variable  to 
flow  problems  has  been  extensively  studied  by  Russian  theoreticians,  particu¬ 
larly  Polubarinova-Kochina  (1962).  Harr  (1962)  gives  an  explanation  of  the 
principles  involved. 


* 
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In  the  analysis  using  the  methods  of  numerical  approximation,  each 
of  the  second  derivatives  in  the  Laplace  equation  (equation  2.1b)  is  expressed 
in  finite  difference  form.  The  total  head  h  at  any  intersection  point,  or 
node,  on  an  imaginary  rectangular  mesh  of  small,  but  finite,  size  is  then  re¬ 
lated  to  the  head  at  the  adjacent  nodes.  An  initial  guess  is  made  of  the  value 
of  the  head  at  each  node  and,  by  using  the  finite  difference  equation,  a  better 
approximation  to  the  correct  value  can  be  made.  The  process  is  repeated  until 
the  difference  in  values  obtained  by  two  successive  approximations  for  all 
nodes  is  less  than  some  tolerable  value.  The  method,  slightly  modified,  is 
sometimes  called  'relaxation',  in  which  an  assessment  is  made  of  the  amount  by 
which  the  guessed  head  at  any  node  is  out-of-balance  (in  terms  of  the  finite 
difference  equation)  with  the  heads  at  adjacent  nodes.  The  guessed  head  is 
then  corrected  by  that  amount,  or  relaxed;  this  correction  throws  the  adjacent 
nodes  out-of-balance,  so  then  they  are  relaxed,  and  the  process  repeated  until 
the  amount  of  unbalance  is  tolerably  small.  Numerical  approximation  and 
relaxation  methods  were  developed  largely  by  R.  V.  Southwell  (see,  for  example, 
Shaw  and  Southwell,  1941,  and  Southwell,  1946). 

(b)  Experimental  methods 

Perhaps  the  method  most  widely  used  by  engineers  for  the  solution  of 
flow  problems  is  the  graphical  construction  of  the  flow  net.  This  method  was 
devised  by  P.  Forchheimer  in  the  early  part  of  the  twentieth  century,  and  con¬ 
sists  of  sketching  the  flow  net,  starting  with  the  known  boundary  conditions, 
proceeding  by  trial-and-error ,  and  recognizing  the  fact  that  equipoten tial  lines 
and  flow  lines  are  orthogonal,  so  that  the  final  picture  consists  of  a  system 
of  'curvilinear  squares'.  This  method  is  described  in  detail  by  A.  Casagrande 


. 


(1937).  Typical  flow-nets  which  are  useful  as  a  guide  have  been  presented 
by  Gregg  (1940). 
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Hydraulic  models  may  be  used  to  determine  the  flow  conditions  in 
or  around  an  earth  structure.  The  structure  is  constructed  to  a  suitable 
scale  between  two  parallel  glass  plates  so  that  the  flow  occurs  parallel  to 
the  plates.  A  head  of  water  is  applied  to  the  upstream  side  of  the  structure 
and  when  flow  has  reached  a  steady-state  condition,  continuous  streams  of  a 
suitable  dye  are  introduced  at  various  points  on  the  upstream  face,  against 
the  side  of  the  glass.  The  traces  of  the  dye  follow  flow  lines,  the  signifi¬ 
cant  coordinates  of  which  may  be  recorded  or  plotted.  Piezometers,  in  the  form 
of  small  bore  tubes,  may  be  inserted  into  the  flow  region  at  specific  points 
and  the  total  head  at  those  points  determined  from  readings  on  manometers  con¬ 
nected  to  each  piezometer.  Knowing  the  flow-lines  and  the  head  at  some  points, 
the  equipotential  lines  may  by  sketched  in.  Taylor  (1948,  FIGURE  9.15)  shows 
an  example  of  a  model  dam.  Zangar  (1953,  FIGURE  38)  shows  photographs  of  flow¬ 
lines  under  a  model  weir  and  under  a  cut-off  wall. 

(c)  Analogue  methods 

The  Laplace  equation,  as  well  as  applying  to  seepage  problems,  rep¬ 
resents  the  condition  governing  potential  distribution  in  several  other  physical 
phenomena.  The  flow  of  electricity  or  heat  through  a  conducting  medium  are 
examples.  The  same  form  of  equation  describes  the  displacement  of  a  loaded 
membrane  perpendicular  to  its  plane;  also,  the  capillary  flow  of  a  liquid  be¬ 
tween  two  closely  spaced  parallel  plates  satisfies  the  condition  expressed  in 
equation  (2.1b).  Some  consideration  and  effort  has,  in  the  past,  been  given 
to  the  use  of  these  phenomena  for  solving  seepage  problems  by  analogy. 
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When  electricity  passes  through  a  two-dimensional  medium  with 
resistivities  Rx,y  in  the  x-  and  y-  directions  respectively,  the  following 
law  is  obeyed: 


where  V  is  the  voltage  at  an  arbitary  point  in  the  plane. 

In  the  analogue,  it  is  necessary  to  obtain  an  electrical  conductor  having 
directional  resistivities  inversely  proportional  to  the  corresponding 
directional  permeabilities.  If  a  rectangular  mesh  with  sides  parallel  to  the 
directions  of  principal  permeabilities  is  imagined  to  be  superimposed  on  the 

prototype,  the  nodes  of  the  corresponding  mesh  in  the  analogue  may  be  formed 

by  the  intersection  of  electrical  resistors  related  to  the  permeabilities. 
Alternatively,  if  the  soil  prototype  is  isotropic,  a  suitable  electrically 
conducting  paper  may  be  used  in  the  analogue.  The  model  is  made  geometrically 
similar  to  the  prototype,  and  uniform  voltages  are  applied  along  each  of  the 

boundaries  representing  faces  of  constant  total  head.  The  voltage  at  any 

point  in  the  paper  analogue  or  any  node  in  the  mesh  analogue,  is  then  measured 
by  use  of  a  probe  and  expressed  as  a  proportion  of  the  total  voltage  difference. 
In  this  way,  a  picture  of  the  potential  distribution  is  built  up  throughout  the 
flow  region.  Cases  for  which  the  electric  analogue  have  been  used  are  shown 
by  Zangar  (1953). 

Although  the  laws  governing  hydraulic  and  thermal  flows  are  analogous, 
solutions  to  hydraulic  flow  problems  by  analogy  with  thermal  models  are  not 
attempted  because  of  difficulties  in  setting  up  the  apparatus  and  in  measur¬ 
ing  the  temperatures  during  the  conduct  of  the  test.  Usually,  the  solution 
of  a  heat  flow  problem  poses  as  large  a  question  as  the  solution  of  a  hydraulic 
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flow  problem. 


If  some  part  of  a  uniformly  stretched  membrane  is  displaced  a 
small  distance  perpendicular  to  its  original  plane  of  location,  the  dis¬ 
placement,  z,  of  any  point  from  that  plane  is  given  by  i 


dzz 
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where  x^and  y-)  are  perpendicular  directions  in  the  original  plane 
of  location. 

A  membrane  having  the  same  shape  as  the  soil  prototype  is  uniformly  stretched 
and  displaced  at  the  edges  by  amounts  proportional  to  the  total  heads  acting 
at  the  corresponding  edges  in  the  prototype.  The  displacement  of  any  point 
on  the  membrane  will  -then  be  proportional  to  the  head  at  the  corresponding 
point  in  the  prototype.  The  method  is  limited  to  cases  of  homogeneous  soils, 
and,  so  far  as  the  writer  is  aware,  it  has  been  little  used  for  seepage  prob¬ 
lems.  For  a  fuller  description  of  the  method,  reference  should  be  made  to 
Walker  (1949),  who  used  a  rubber  membrane  analogue  to  investigate  the  motion 


of  charged  particles  in  an  electric  field. 


In  the  capillary  flow  analogue  a  viscous  fluid  is  allowed  to  pass 
between  two  closely  spaced  parallel  plates  which  have  the  same  shape  and  boun¬ 
dary  conditions  as  the  soil  profile.  The  spacing  of  the  plates  is  determined 
in  conjunction  with  the  viscosity  of  the  fluid  so  that  the  Reynold*s  number  is 
subcritical  and  laminar  flow  exists.  Polubarinova-Kochina  (1962,  pp. 465-477) 
describes  the  method  and  shows  how  nan-homogeneous  soils  may  be  simulated  in 
the  analogue. 


2.3  Appraisal  of  existing  methods 


Because  seepage  problems  occupy  the  interests  of  several  groups  of 
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persons,  such  as  soil  engineers,  petroleum  engineers,  mathematicians  and 
physicists,  there  exists  an  abundance  of  material  in  the  engineering 
literature  dealing  with  methods  of  analysis.  A  few  references  are  made 
to  specific  texts  in  this  section;  for  further  material,  the  bibliographies 
in  the  works  of  Scheidegger  (1957)  and  Harr  (1962)  should  be  consulted. 

The  method  of  analysis  involving  the  graphical  construction  of 
flow  nets  has  several  advantages;  no  equipment  or  specialised  mathematical 
knowledge  is  required,  it  is  simple  in  operation,  and  as  the  user  progresses 
through  each  trial-and-error  stage  of  the  operation,  he  develops  an  intuitive 
understanding  of  the  conditions  which  are  most  likely  to  affect  the  flow  net. 
The  method  has  received  severe  criticism  from  Leliavsky  (1955),  who  believes 
that  the  basic  criterion  of  forming  'curvilinear  squares'  was  not  rigid 
enough;  he  suggests  that  better  'squares’  could  be  formed  by  constructing  an 
array  of  circles  inside  the  flow  region  so  that  adjacent  circles  have  common 
tangent  points  (as  shown  for  the  flow  region  ABCD  in  FIGURE  2.1).  However, 
since  a  poorly  drawn  flow  net  usually  produces  no  more  than  a  10  per  cent 
error  in  the  calculation  of  the  seepage  quantity,  Leliavsky's  criticism  hard¬ 
ly  seems  warranted.  Two  main  difficulties  are  experienced  if  the  flow  region 
consists  of  a  non-homoqeneous  soil;  first,  although  a  relationship  exists  for 
the  change  of  direction  of  flow  lines  at  their  intersection  with  a  common 
boundary  between  two  zones,  the  expression  is  awkward  to  use  in  the  graphical 
method;  secondly,  if  the  zones  have  different  degrees  of  anisotropy  the  stan¬ 
dard  scale  transformation  (Taylor,  1948,  chapter  9)  will  distort  some  zones 
more  than  others  and  corresponding  points  on  zone  boundaries  might  not  remain 
coincident.  Barron  (1948)  has  presented  a  technique  which  reduces  these 
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(a)  Construction  of  circles  with  common  tangent  points. 


(b)  Insertion  of  flow  lines  and  equipotential  lines. 


(c)  Comparison  with  constructed  flow  net  (full  lines) 
with  sketched  flow  net  (broken  lines). 


Circle  method  of  flow  net  construction 

(after  Leliavsky,  1955) 


FIGURE  2.1. 
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difficulties,  but  the  quantity  of  work  required  for  a  solution  of  the  case 
of  a  non-homogeneous,  anisotropic  soil  by  graphical  construction  still  re¬ 
mains  large. 

Hydraulic  models  are  an  ideal  means  by  which  the  flow  net  in  the 
prototype  can  be  visualized.  For  flow  regions  having  a  phreatic  line  as  one 
boundary,  the  use  of  models  enables  this  line  to  be  determined  directly;  how¬ 
ever,  as  Taylor  (1940)  points  out,  capillarity  effects  may  detract  from  the 
accuracy  of  the  determination  of  this  line.  In  order  to  reduce  the  capillar¬ 
ity  effects  in  the  model,  a  coarse  sand  is  often  used.  Even  so,  as  Kerr  (1959) 
concludes  from  some  earth  dam  model  tests  using  sand,  the  effect  of  the  capillary 
fringe  may  not  be  negligible.  As  well  as  the  difficulty  in  locating  the  capil¬ 
lary  fringe  accurately,  capillarity  effects  have  a  further  manifestation; 
there  exists  a  variation  in  degree  of  saturation  of  the  soil  within  the  fringe. 
Because  the  permeability  to  water  of  a  soil  is  dependent  on  its  degree  of  satur¬ 
ation,  permeability  should  strictly  be  treated  as  a  variable  in  flow  analyses. 
Liakopoulos  (1964)  has  conducted  tests  on  a  soil  inside  a  tube  open  at  the  ends 
and  through  which  the  flow  was  one-dimensional;  from  these  tests,  the  relation¬ 
ships  between  degree  of  saturation,  permeability,  and  position  inside  the  capil¬ 
lary  fringe  were  determined.  The  relationships  appertain  to  any  seepage  con¬ 
dition  satisfying  Darcy's  law,  but  the  difficulty  lies  in  applying  them  to  two- 
dimensional  flow  conditions.  Because  of  the  doubt  which  exists  on  the  true 
effect  of  capillarity,  it  is  usual  to  assume  that  the  phreatic  line  lies  along 
the  upper  limit  of  saturation. 

The  parallel  plate  analogue  permits  direct  determination  of  the  phreatic 
line;  similation  of  the  capillary  fringe  in  the  prototype  may  be  made  in  the 
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analogue  by  correct  spacing  of  the  plates,  as  described  by  Polubarinova- 
Kochina  (1962).  Disadvantages  of  the  parallel  plate  analogue  are  that 
measurement  of  the  potential  head  at  any  point  is  practically  unobtainable, 
and  that  if  the  plates  are  not  perfectly  parallel,  large  errors  will  result. 

The  electric  analogue  method  has  more  common  use  than  other  methods, 
since  it  has  several  advantages;  as  well  as  being  simple  and  convenient,  steady- 
state  flow  occurs  immediately,  so  that,  once  the  equipment  is  set  up,  the  pro¬ 
cedure  is  usually  completed  quickly.  One  disadvantage  is  that  phreatic  lines 
must  be  located  by  trial-and-error .  The  existence  of  the  capillary  fringe 
has  been  taken  into  account  in  tests  using  the  electric  analogue  by 
Vreedenburgh  (1936)  and  Chapman  (i960). 

The  main  problem  with  the  complex  variable  method  of  flow  analysis, 
as  Harr  and  Deen  (1961)  point  out,  lies  in  finding  an  equation  which  will  per¬ 
form  the  desired  transformation.  Any  polygonal  region  in  one  plane  may  be 
mapped  on  to  the  upper  half  of  an  auxiliary  plane  by  means  of  the  Schwarz- 
Christoffel  transformation  (Milne-Thomson ,  I960).  Although  the  transformation 
then  breaks  down  into  a  step-by-step  procedure,  it  might  be  complicated  for 
any  case  other  than  the  simplest.  Because  of  this  difficulty,  the  method  is 
of  limited  value  to  soil  engineers,  even  though  published  solutions  by  this 
method  to  certain  practical  cases  ( Polubarinova-Kochina ,  1962)  and  lists  of 
conformal  transformations  (for  example,  Kober,  1952)  are  available.  For  cases 
in  which  the  boundary  conditions  are  simply  described,  the  method  is  a  mathe¬ 
matical  exercise  of  reasonable  length.  A  classical  example  is  the  case  attempted 
by  Pavlovsky  (1936).  He  develops  the  solution  to  flow  under  an  impervious 
structure  founded  on  the  surface  of  a  homogeneous  permeable  soil  of  infinite 
depth.  From  his  solution,  an  expression  for  the  distribution  of  uplift  pressure 
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across  the  base  of  the  structure  can  immediately  be  derived.  Similarly, 
the  complex  variable  theory  may  advantageously  be  applied  in  the  case  of 
flow  around  sheet-piling  penetrating  some  distance  into  a  homogeneous  perm¬ 
eable  soil  of  infinite  depth,  to  determine  an  expression  for  the  upward 
hydraulic  gradient  near  the  exit.  It  might  be  argued  that  these  conditions 
are  idealistic  but  they  form  a  good  basis  for  making  a  prediction  of  some 
of  the  conditions  in  more  realistic  problems. 

At  the  present  stage  of  engineering  science,  it  would  be  an  ill- 
advised  procedure  to  attempt  a  solution  to  the  case  of  flow  through  a  non- 
homogeneous  anisotropic  soil  by  use  of  complex  variable  theory.  Further¬ 
more,  as  already  explained,  such  a  problem  would  be  cumbersome  to  solve  by 
the  graphical  method  and  difficulties  of  a  practical  nature  would  arise  if  a 
model  or  analogue  were  used.  It  is  possible  to  solve  the  problem  using  the 
conventional  methods  of  numerical  approximation  but,  even  so,  difficulties 
are  encountered.  These  difficulties  are  explained  in  the  next  section,  to¬ 
gether  with  the  writer's  proposal  for  overcoming  them. 

2.4  Numerical  analysis  of  flow  problems 

For  the  sake  of  brevity,  the  following  description  of  numerical 
approximate  methods  considers  flow  in  two  dimensions  only,  through  an  iso¬ 
tropic  soil.  The  order  of  complexity  would  not  be  enlarged  if  the  third  dimen¬ 
sion  and  anisotropy  were  included. 

The  Laplace  equation  (equation  2.1b)  can  be  stated  in  the  following 
finite  difference  form: 

h,  +  h*  +  h3  +  h*  -  ^ho  “  O,  (2.2a) 


.O  -  0rrt*  ■*  gH  *  ■*  ,r\ 
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where  hD  is  the  total  head  at  a  node  in  an  imaginary  square 

mesh,  with  small  sides,  superimposed  on  the  flow  region, 
and  are  the  total  heads  at  each  of  the  four  adjacent  nodes. 


The  derivation  of  equation  (2.2a)  may  be  found  in  most  textbooks 
on  rela>atiam  methods  (for  example,  Shaw,  1953,  chapter  2).  Pictorial  rep¬ 
resentation  of  the  mesh  is  shown  in  FIGURE  2.2(a).  By  ensuring  that  equation 
(2.2a)  is  satisfied  at  all  nodes  in  the  mesh,  an  approximation  is  obtained 
to  the  distribution  of  total  head  within  the  flow  region.  Generally,  the 
boundaries  of  the  flow  reqion  and  of  the  individual  zones  do  not  lie  along 
a  mesh  side  (as  shown  in  FIGURE  2.2(b)),  and  either  an  approximate  boundary 
may  be  drawn,  following  the  true  boundary  as  closely  as  possible  but  keeping 
strictly  to  mesh  lines,  or  a  modification  to  equation  (2.2 a)  is  required  for 
those  nodes  adjacent  to  the  boundary.  The  modified  finite  difference  equation 
for  nodes  near  boundaries  is  (Shaw,  1953,  chapter  5): 


Ah,  4-  Bh2  4  Ch3  4  Dh4  -  (E4F)h0  »  O,  (2.2b) 

where 
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al,a2  are  the  lengths  shown  in  FIGURE  2.2(b), 

the  locations  of  and  T14  are  as  shown  in  FIGURE  2.2(b) 

that  for  each  node  less  than  ohe  mesh  length  from  the  boundary 

coefficients  A,....,F  in  equation  (2.2b)  are  generally  dif- 

on  the  horizontal  and  vertical  distances  from  the  node  to 


The  techniques  used  in  numerical  approximation  have  already  been 
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(a ) 


FIGURE  2.2. 


Examples  of  the  square  relaxation  mesh. 


20 


mentioned  (  5ECTI0N  2.2  ) .  Inspection  of  equations  (2.2a)  and  (2.2b)  indicates 
that  nodes  not  adjacent  to  boundaries  are  no  cause  for  difficulty,  but  nodes 
adjacent  to  boundaries  require  individual  treatment;  thus  if  the  boundaries 
are  not  approximated  along  the  mesh  sides,  the  method  becomes  appreciably  more 
laborious  if  being  performed  by  hand,  or,  if  an  electronic  computer  is  being 
used,  so  much  more  programming  is  required. 

Gibson  and  Lumb  (1953)  appear  to  have  reduced  the  difficulty  of 
nodes  not  lying  on  the  boundaries;  in  an  analysis  of  an  earth  dam,  the  shape 
was  transformed  so  that  the  base  angles  became  45°;  a  square  mesh,  with  nodes 
lying  on  the  inclined  faces,  was  used.  5cott  (1963,  chapter  4)  has  suggested 
that,  for  the  case  of  a  triangular  dam,  a  mash  consisting  of  equilateral  tri¬ 
angles  might  advantageously  be  used.  By  transforming  the  vertical  scale  of 
the  dam,  the  faces  could  be  arranged  to  lie  at  60°  to  the  horizontal.  The  en¬ 
suing  computations  would  then  be  less  laborious. 

The  writer's  proposal  for  dealing  with  the  problem  of  flow  through 
a  non-homogeneous ,  anisotropic  soil  was  to  take  Scott's  suggestion  of  a  tri¬ 
angular  mesh  even  further.  Having  approximated  any  curved  boundaries  that 
are  present  in  the  flow  region  to  a  series  of  straight  lines,  the  section  is 
divided  into  triangles  so  that  each  zone  boundary  lies  along  a  side  of  one  of 
the  triangles.  The  flow  region  can  then  be  regarded  as  looking  like  a  jig-saw 
puzzle  in  which  the  pieces  consist  of  a  variety  of  dissimilar  triangles.  Each 
triangle  is  then  sub-divided  into  smaller,  similar  triangles  by  drawing  equally 
spaced  lines  parallel  to  each  of  the  sides.  The  procedure  is  illustrated  in 
FIGURE  2.3.  By  using  this  procedure  for  forming  a  mesh,  the  subsequent  num¬ 
erical  work  is  systematized,  since  the  positions  nf  nodes  always  fall  into  one 
of"  three  categories,  as  follows!  (a)  inside  a  large  triangle  and  not  less  than 
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(a)  Flow  region  with  arbitary  boundaries  AB  and  CD 
inserted  upstream  and  downstream. 


(b)  Flow  region  divided  into  dissimilar  triangles. 


(c)  Mesh  configuration;  for  clarity,  each  triangle  is 
shown  sub-divided  into  only  nine 
smaller  triangles;  actually  there 
may  be  any  number  of  sub-divisions. 


FIGURE  2.3. 


Procedure  for  forming  mesh. 


. 
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one  mesh  length  distant  from  any  of  its  edges,  (b)  on  an  edge,  but  not  at 
an  apex,  of  a  large  triangle,  or  (c)  at  an  apex  of  a  large  triangle. 

As  is  indicated  in  the  next  chapter,  the  finite  difference 
equations  associated  with  the  proposed  type  of  mesh  are  more  complicated 
than  the  equation  used  with  a  square  mesh  (equation  2.2a),  It  is  this 
additional  complexity  that  makes  a  hand  solution  extremely  laborious,  but 
causes  little  extra  difficulty  when  a  computer  is  available. 


CHAPTER  III 


DERIVATION  OF  FINITE  DIFFERENCE  EQUATIONS 

3,1  Introduction 

Most  textbooks  which  give  the  derivation  of  the  'basic  equation 
of  flow  do  so  by  considering  an  elemental  rectangle,  or  rectangular  prism 
in  the  case  of  three-dimensional  flow,  having  sides  parallel  to  the  direc¬ 
tions  of  principal  permeability.  The  derivation  uses  Darcy's  law  to  obtain 
an  expression  for  the  flow  across  each  face  of  the  rectangle  in  terms  of 
the  hydraulic  gradient;  the  algebraic  sum  of  these  expressions  is  made  equal 
to  zero  in  order  to  satisfy  the  continuity  condition.  As  explained  in  SECTION 
2.4,  a  triangular  mesh  is  adopted  herein,  so  that  difficulties  at  zone  bound¬ 
aries  are  avoided;  this  mesh  configuration  has  three  axes,  each  parallel  to 
one  of  the  sides  of  the  triangle.  The  angle  between  each  axis  and  some  datum 
direction  has  any  general  value,  making  it  advantageous  to  use  a  circle  in¬ 
stead  of  a  rectangle  for  the  purpose  of  deriving  a  continuity  equation.  If 
this  circle  is  centred  on  the  intersection  of  the  axes,  each  axis  becomes  a 
radius  and  some  symmetry  is  retained,  making  the  algebra  required  in  the 
derivation  somewhat  easier  than  if  a  rectangle  had  been  used.  Apart  from 
Darcy's  law  and  the  continuity  condition,  no  additional  physical  concepts  are 
introduced  in  the  derivation  of  the  continuity  equation  in  terms  of  three  axes 
in  a  plane. 

The  finite  difference  equations  are  a  special  form  of  the  continuity 
equation;  they  are  partly  obtained  through  approximate  expressions  for  the 
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derivatives  of  head  with  respect  to  distance.  The  equations  always  appear 
herein  in  a  form  relating  the  head  at  an  arbitary  node  D  to  the  adjacent 
nodes  on  all  the  mesh  axes  through  0. 

The  following  sections  describe  the  principal  steps  in  the  deriva¬ 
tion  of  the  continuity  equation  and  the  finite  difference  equations.  Greater 
detail  of  the  derivations  is  presented  in  appendix  A  (continuity  equation)  and 
appendix  B  (finite  difference  equations). 

3.2  Continuity  enuation 

The  finite  difference  equations  for  each  category  of  node  cannot  be 
obtained  until  the  continuity  equation  in  terms  of  three  axes  in  a  plane  has 
been  derived. 

Consider  the  elemental  circle  centred  on  an  arbitary  point  0  in  an 
anisotropic  soil  having  principal  permeabilities  kx,y  in  the  x-  and  y-  direc¬ 
tions  respectively  (FIGURE  3.1).  By  consideration  of  rates  of  change  of  super¬ 
ficial  velocity  (hereinafter  called  simply  'velocity'),  an  expression  is  ob¬ 
tained  for  the  outward  flow  q  across  the  elemental  arc  PP' ,  in  terms  of  the 
velocity  components  v  ,  at  0  in  the  x-  and  y-  directions  respectively,  as 

x  y 

follows^  : 

<\  =  {  Vx  cos  4-  Vy  Sin  * }  .  As  . 

+  {  (vx)3  COS*  4-  (vy)5  si  (As)*  Sac 

2 

Equation  (3.l)  can  be  altered  to  equation  (3.3)  by  use  of  Darcy's 

law. 

]See  footnote, appendix  A,  page  A2. 

^"Equations  in  this  chapter  are  numbered  to  agree  with  a  ppendices  A  and  B,  in 
which  they  are  numbered  consecutively. 


'  ' 


FIGURE  3.1. 


Elemental  circle. 


^  {  Vx  CoS  <*  -f  Vy  sin  }  .  As.  <£*  I 

~  {  (h)xs  CoSot  4-  ky  (k)ys  Sin  *  j- .(As)*  <£*  J 

The  continuity  equation  is  conveniently  expressed  in  terms  of  the 

three  axes  Ou,  Ov  and  Ow.  The  terms  (h)  and  (h)  in  equation  (3.3)  can  be 

x  s  vs 

expressed  in  terms  of  the  second  derivatives  of  h  along  any  two  of  the  Ou,  Ov 

and  Ow  axes.  The  third  axis  is  used  to  eliminate  the  term  comprising  the 

mixed  second  derivative  of  h  along  the  first  two  axes.  For  example,  if  Ov 

is  chosen  as  the  third  axis,  (h)  and  (h)  are  expressed  as  functions  of 

(h)  ,  (h)ww  and  (h)uw*  Another  expression  relating  (h)w,  ( h ) uu,  (h)ww  and 

(h)uw  is  used  to  eliminate  (h)uw,  thus  leaving  (h)xs  and  (h)  s  expressed  in 

terms  of  (h)  ,  (h)  and  (h)  .  As  would  be  expected,  these  expressions  are 

uu'  w  ww 

symmetrical  in  u,  v  and  w  and  may  be  written  in  the  convenient  forms  of  equations 
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(3.9c)  and  (3.9d). 


(h)*S  =  y1  (^MM  '  S'0  ~Sir>  («-*?>+  sinolP-  5'"6*-«Nl)  ( 3,9c , 

M=^w  2  S,Vl  C“n-'m)'5Io(o(M-<p) 


and 


(V)  =  V  M  .  CoS<tfN.  Sin  (*-<*p)  +  CoSo<p  Slo(o(-o(N) 

V  "P  /_i  '  '  wm  - o~~7' — 7 - ^  N  .  > - ^ — 

where  M,  N  and  P  are,  respectively,  u,  v,  w, 

and  v,  w,  u, 
and  w,  u,  v. 


( 3 . 9d  ) 


Equations  (3.9c)  and  (3.9d)  may  be  substituted  in  equation  (3.3), 
which  is  then  rearranged  to  give  equation  (3.10)* 


(J,*  jvrGOSO(  +  Vy  sin  As.  Su 

JV  /l]  j  (kx+  ki4)sin(c(N-»-c<p)ce6ctSinc<  -2kx$IO  attain  -  2k.a6o«MC^6Wr5intB<7l|lA^. 

(H'  2  »n  (*„  -  «tM)  5m  (ctM-  otP)  jJ  y 


(3.10) 


Continuity  requires  that  the  net  outflow  from  the  circle  is  zero 
when  0  is  a  point  inside  a  homogeneous  soil.  This  condition  leads  to  equation 
(3.11). 


L 


Th')  .  k*  Sin  o<N  Sin  ocp  +  ky  CoS  CoS  «(,  s  q  (3.11) 

t  /MM  - 71 — 7 - c - : - - ; - : - 

5m  ~  «*w  )  •  6»n  (o(|4  —  otp) 


M=  u,v,W 


Equation  (3.11)  is  the  continuity  equation  of  flow,  in  terms  of  three 


axes  in  a  plane. 


3.3  Finite  difference  equation:  point  inside  a  homogeneous  soil 

For  steady  state  flow  conditions,  the  continuity  equation  is  satis¬ 
fied  at  all  points  in  the  flow  region.  When  a  mesh  is  superimposed  on  the  flow 
region ,  an  approximation  to  the  steady— state  flow  conditions  is  obtained  if  the 
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FIGURE  3.2  Part  of  a  triangular  mesh. 

continuity  equation  is  satisfied  at  all  nodes.  The  difference  between 
the  approximation  and  the  true  conditions  is  small  if  a  fine  mesh  is  chosen. 

The  concept  of  a  mesh  is  useful  as  a  means  for  expressing  the  sec¬ 
ond  derivatives  of  h  appearing  in  the  continuity  equation  (equation  3.11)  in 
a  form  which  is  easily  soluble.  For  example,  at  point  0  in  FIGURE  3.2,  the 
second  derivative  of  h  with  respect  to  distance  measured  along  the  Ou  axis 
can  be  shown  to  be  approximately  equal  to  the  right-hand  side  of  equation 
(3.14a). 

+  hu  "  2h0  i  , 

V'l/uu  - 7~ - nj -  (3.14a) 

(A»«r 

where  h^,  h",  hQ  are  the  total  heads  at  U',  UM,  0  respectively, 
and  A|U  is  the  length  QU*. 

The  same  form  of  finite  difference  expression  can  be  obtained  for  (h)  and 
(h)ww  at  the  point  0. 


Referring  to  FIGURE  3.2,  it  may  be  observed  that,  by  the  sine  rule 


for  triangles, 
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A,u  _  A,  v  _ 

Sin  («w— olv)  jin  (oiw  -ct  u) 

thus  sines  of  angles  can  be  substituted 
denominators  of  the  finite  difference  expressions 


_ A  |  W _  • 

Sm(*v-c<u) 

for  the  lengths  in  the 
such  as  equation  (3.14a). 


The  finite  difference  expressions  may  then  be  substituted  in  equation 
(3.11),  giving 

T  (C+ln|!,-2h(,).  5'"  “n  sin  rfP  +  kj,  o>s  *N  cos«,  _  q  ,  (3>17) 

M=u,v,w  Sin  C«n-«p) 

which  can  be  written  as 

+  Cy.  hy  +  Cw*hw  +■  Cu*^U  +•  Cyhy  +  Cw  “2  (Cu  +Cy  +■  Cw)h0  =  O,  (3.10) 


where 


c  =  kx  sin  (<*i -olft)  sin  (oCyl/ +  ky  CoS  ~^a)  (3.19a) 

U  '  /  I  I  \  J 

Sll^  ^oCy  "*  Ctyy^  J 

Cy  -  ^<3C  Sin  (o^w -opsin  kg  oaS  (^w-o^g)  OaS  (o<u -eta)  >  (3  19bj 

sin  (oijy  ~  o(  u) 

Cv  =  k*'  6<v-o<o)  +  0X3  (o(l  -0(a)  Cos  (oCy  -  o(a)  %  (3  19c) 

Stn  (o<w  -  o<a) 

and  <*' ,  oc'  ,  o<* »  are  the  angles  measured  in  a  counter-clockwise 
u  v  w 

sense  from  the  datum  (horizontal)  direction  to  the  axes  Ou, 
Ov,  Ow  respectively. 


The  coefficients  C^,  C^,  and  may  be  computed  for  any  given  data, 


w 


and  hQ  may  be  written  in  terms  of  the  heads  at  surrounding  nodes  by  using 


equation  (3.10), 
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FIGURE  3.3.  Point  on  a  boundary. 

3.4  Finite  difference  equation:  point  on  a  boundary  between  two  zones 

All  nodes  on  common  boundaries  of  two  zones  are  common  to  the 
meshes  of  both  zones,  because  the  sides  of  each  triangle  formed  by  dividing 
the  polygonal  zones  are  sub-divided  equally  into  the  same  number  of  segments. 
A  boundary  is  any  one  of  the  Ou,  Ov  or  Ow  axes  in  terms  of  the  mesh  config¬ 
uration  of  one  of  the  zones;  it  is  also  any  one  of  the  Ou,  Ov  or  Ow  axes  in 
terms  of  the  other  zone.  FIGURE  3.3  shows  the  Ow  axes  of  each  zone  lying 
along  the  boundary  (depicted  by  the  double  line).  The  derivation  appearing 
in  this  section  applies  to  the  case  as  shown. 

The  axes  which  do  not  lie  on  the  boundary  (i.e.,  the  Ou  and  Ov  axes 
in  FIGURE  3.3)  do  not  generally  have  the  same  alignment  either  side  of  the 
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boundary,  and  therefore  a  finite  difference  expression  having  the  same  form 
as  equation  (3.14a)  cannot  be  stated  directly  for  these  axes.  In  order  that 
the  concept  of  finite  differences  may  still  be  used,  each  zone  is  considered 
separately,  and  an  image  of  the  mesh  is  produced  on  the  opposite  side  of  the 
boundary.  In  FIGURE  3.4,  one  of  the  zones  shown  in  FIGURE  3.3  is  redrawn, 
with  its  image.  The  subscript  1  denotes  the  zone  to  which  the  notation  refers. 
Ficticicus  heads  are  introduced  at  the  image  points  U^  and  so  that  finite 
difference  expressions  similar  in  form  to  equation  (3.14a)  can  be  stated. 

f  f 

The  ficticious  heads  at  U^  and  and  the  corresponding  ficticious 
heads  for  zone  2  comprise  four  unknowns.  Five  simultaneous  equations,  each 
defining  some  independent  physical  aspect  of  the  flow  conditions  in  the 
vicinity  of  point  0,  are  therefore  required;  four  equations  are  used  in  the 
elimination  of  the  four  unknowns,  and  the  remaining  one  becomes  the  finite 
difference  equation.  The  simultaneous  equations  are  obtained  from  the  fol¬ 
lowing  conditions: 

(a)  the  continuity  of  flow  within  the  elemental  semi-circle  GH  centred 
on  0,  in  zone  1,  (FIGURE  3.4), 

(b)  the  continuity  of  flow  within  a  similar  semi-circle  in  zone  2, 

(c)  the  fact  that  only  two  axes  are  required  to  locate  a  point  in  a 

plane  and  therefore  a  distance  measured  along  any  one  of  the  0u^ ,  0v1  or  0w^ 
axes  may  be  regarded  as  a  function  of  the  other  two  axes, 

(d)  the  same  as  (c),  but  using  the  0u2»  0v2  and  0w2  axes, 

(e)  the  flow  into  zone  1  across  GH  is  equal  to  the  flow  out  of  zone  2 

across  GH. 


The  continuity  of  flow  within  the  semi-circle  GH  centred  on  0  is 


Uj  axis 
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FIGURE  3.4  Image  of  zone  1. 


satisfied  if  the  outward  flow  across  the  curved  portion  of  the  perimeter 
is  equal  to  the  flow  Q2  into  zone  1  across  GH.  Integration  of  equation  (3.10) 
between  limits  of  (t*w^  -  IT )  and  otw^  shows  that 


V^,  Sinc*w,  -  V4|  c©6o<w| 
A 


}■  A« 


(3.20) 


Q,  =  2{ 

4  TT  y*  (l\  kxi  Si'rt  o(  Wl  Sin  c/p,  +  ky,  CoSc<?, 

2  Zj  MM|*  Sin  (Wmi-^Pi)  ^ 

where  M,  IM  and  P  have  the  same  meaning  as  for  equation  (3.10), 
The  flow  Q2  across  GH  into  zone  1  may  be  shown  to  be 


S  in 


CoS  oi  \fj | 


(3.22) 
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Since  «  Q2,  it  follows  that 


2,  ‘ 


|  (3. 


24) 


rMM  ._‘X1  Sl0c(N|  °<Pi  +  CcS  <*ni  *>t  O.  (3.23) 

sin  (c^m,-o<pO 

Equation  (3.23)  is  written  in  finite  difference  form  as  follows: 

Cu  I*  ^U|  4-  Cy, •  hy  j  +■  Cw,.  I  •  h  Ul  4*  C  y  1 .  +■  CW|.  ^ W2. 

”"2  CCu»  +  Cv«  ■+■  C'wi)  ho  *  0> 

where  h^  ,  h^  ,  h^  ,  hV,  are  the  total  heads  at  Uj ,  Vj  ,  Wj ,  Wj, 
respectively, 

f  f  f  f 

hu1,  hv^  are  the  ficticious  heads  at  Uj  ,  V-j  respectively, 

and  ,  Cv1 ,  are  the  coefficients  appertaining  to  zone  1 

having  the  same  expansions  as  written  in  equations  (3.19a), 

(3.19b)  and  (3.19c). 


Similarly,  for  zone  2,  the  finite  difference  equation  is 

Cu-2.  hu*2  +■  ^MT.'  2  +  •  k^c„.hf1  + 

Cwz-  h  wi 


"2  (Cu-2.  +  +  CW2,)ln0  ^  O 


.K 


25) 


If  the  Ou^  and  Ow^  axes  are  chosen  for  defining  the  plane  of  flow 
in  zone  1,  (h)v-j  can  be  expressed  in  terms  of  (h)ui  and  (h)w^  as  follows: 

(h)yi  =  (h)u,.(u,)v,  +  ( h)wi- (wi)vi . 


(3.26) 


The  derivatives  of  h  may  be  written  in  finite  difference  form;  for  example, 


for  (h)u1> 


OOul  =  ,  A£-  . 

2(A,u.) 


(3.28a) 


Expressions  such  as  the  right-hand  side  of  equation  (3.28a)  may  be  substituted 
in  equation  (3.26),  leading  to  the  following  relationship: 

h^<  —  hj,  +*  hw,  -  h 4-  hyj  -  ^  .  (3.30) 
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Equation  (3,31)  can  similarly  be  obtained  for  zone  2, 

H “  h  V2.  W2.  “  +  h  va  *“  ^Wi  ^  O  ,  (3.31) 

An  expression  for  the  flow  into  zone  1  across  GH  has  already  been 
obtained,  (equation  3.22).  by  Darcy’s  law, 

Vjjj  =  •  (h)ael  and  I  “  “  kyl  •  (h)yl  y 

these  expressions  may  be  substituted  in  equation  ( 3 ,22 ),  giving 

Q*  =  ~£{kxi»(^)x.iSinc<W|  “  Cos  U  ,  /^S  .  ( 3.36  ) 

The  terms  (h)x^  and  (h)y^  may  be  replaced  by  expressions  containing  (h)u^ 
and  (h)w^,  leading  to  equation  (3.38), 

Qfc  s  {(Kin  +  Cv\)  +  (hi,  -  hw^).  CVj|.{As/A,W,|,  (3.30) 

where  and  Cv^  have  the  same  meaning  as  for  equation  (3.24). 

A  further  expression  for  Q2  can  be  developed  by  considering  zone  2,  i.e.t 

Q2  =  U2  “  ^  U2  )  •  +  C'V'z)  +  0^w*z.  "  ^As/A,W2j'»  (3.39) 

where  Cu2  and  Cv2  have  the  same  meaning  as  for  equation  (3.25). 
Since  wi  =A^W2»  ^  follows  that 

(Cui  +  Cvi)*(hui~hul)  +  (CuZ  +■  C'vz)*  (^Ul  -  h  Ui) 

"*■  C^vi  ^\jz)  •  ( l^wi  “  ^2)  ~ 


j-  (3.40) 


The  term®  h^„  .  h^,  ,  h^->  and  h^_  are  eliminated  from  equations 

ul  *  vl  »  u2  v2  H 

(3.24),  (3.25),  (3.30),  (3.31)  and  (3.40),  leaving  equation  (3.41)  as  the 
equation  for  hQ  in  terms  of  the  heads  at  surrounding  nodes. 


C0,  *^Ul  l^Y2.  +■  .  (hvY("h  ^*W2.) 

“  (Cu  I  C”V»  ■+*  Cw,  4"  C  4-  Cv2  +  CW2.)  -  O, 


(3.41) 
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where  the  Ow  axes  of  both  zones  lie  along  the  boundary;  in  general, 
when  other  axes  lie  along  the  boundary  the  subscripts  u,  v 
and  w  require  interchanging  appropriately. 

Equation  (3.41)  may  be  interpreted  as  follows:  multiply  the  heads 
at  adjacent  nodes  on  the  boundary  by  half  of  the  sum  of  the  coefficients 
(calculated  from  equations  3.19a,  b,  c)  relevant  to  the  axes  on  the  boundary, 
and  add  this  amount  to  the  sum  of  the  heads  at  the  four  other  adjacent  nodes, 
each  multiplied  by  the  appropriate  coefficient;  the  result  is  equal  to  the  sum 
of  all  the  coefficients  multiplied  by  the  head  at  point  0. 

The  finite  difference  equation  for  a  node  on  an  impermeable  bound¬ 
ary  of  the  flow  region  is  easily  obtained  from  equation  (3.41).  An  imperm¬ 
eable  boundary  may  be  regarded  as  a  boundary  between  two  zones,  one  of  which 
has  zero  permeability.  If  the  zones  are  denoted  by  1  and  2,  and  if  zone  2  is 
impermeable,  inspection  of  equations  (3.19a,  b,  c)  shows  that 

^U2  —  C Y2.  ~  ^  f 

in  equation  (3.41),  Thus  omitting  the  subscript  1,  it  follows  from 
equation  (3.41)  that,  for  a  node  on  an  impermeable  boundary, 

CUX  +  Cv.h;  +  Cw .  hw  i±M  -  (Cu  +  cv  +  Cw).ho*o,  (3.42) 

where  ,  h(J,  represent  the  heads  at  nodes  on  the  boundary  either  side 
of  the  node  at  0. 

3.5  Finite  difference  equation:  point  at  an  apex  of  several  triangles 

In  setting  up  the  mesh  for  a  non-homogeneous  soil,  each  polygonal 
zone  is  divided  into  triangles  by  joining  certain  pairs  of  its  apices.  Since 
the  polygons  are  irregular  in  shape  and  the  number  of  sides  of  each  polygon 
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FIGURE  3.5.  The  apex  of  several  triangles. 

is  not  a  constant,  each  apex  of  any  resulting  triangle  may  be  common  to 
several  triangles,  as  shown  in  FIGURE  3.5. 

Despite  considerable  effort,  a  rigorous  proof  of  the  finite  dif¬ 
ference  equation  for  the  head  at  an  apex  of  several  zones  was  not  found.  The 
most  promising  approach  used  by  the  writer  to  establish  such  an  equation  was 
similar  to  the  method  used  for  obtaining  the  finite  difference  equation  for 
a  point  on  the  boundary  between  two  zones.  An  elemental  circle  was  drawn 
with  its  centre  at  0  and  equations  for  the  continuity  condition  for  the  sec¬ 
tor  of  the  circle  lying  in  each  triangle  were  determined  in  a  finite  differ¬ 
ence  form,  using  ficticious  heads  at  image  points.  Also  determined  were  the 
flow  conditions  across  each  boundary,  and  equations  based  on  the  fact  that  a 
measurement  along  one  axis  of  any  triangle  is  a  function  of  the  other  two 
axes  of  that  triangle.  These  conditions  were  insufficient  to  yield  an  equa¬ 
tion  of  the  desired  form,  and  the  writer's  attempts  to  establish  such  an  equa¬ 


tion  by  rigorous  methods  were  abandoned  in  favour  of  a  more  heuristic  approach 
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Inspection  of  the  finite  difference  equations  for  points  inside 
a  homogeneous  soil  and  on  boundaries  (equations  3.18  and  3.41)  shows  that,  for 
these  cases,  the  value  of  h0  is  a  weighted  average  of  the  heads  at  every  ad¬ 
jacent  node.  Except  when  0  is  joined  to  the  adjacent  node  by  a  boundary,  the 
head  at  the  node  is  multiplied  by  the  coefficient  C  having  the  subscript  of 
the  same  letter  u,  v  or  w  as  the  axes  through  0  on  which  the  node  lies.  If 
0  and  the  adjacent  node  are  joined  by  a  boundary,  the  head  at  the  adjacent 
node  is  multiplied  by  half  of  the  sum  of  the  coefficients  C  relevant  to  the 
axes  of  both  zones  which  lie  along  the  boundary.  Thus,  it  appears  reason¬ 
able  to  infer  that  the  head  at  0,  an  apex  of  several  triangles,  is  the 
weighted  average  of  the  heads  at  adjacent  nodes.  Fox  example,  in  FIGURE  3.5, 


Cv8  |qu,+  ^U^  +  Cy,  ^4  CU3  +  CV2  Cu4+  Cy3  Cu54CV4  h,(, 

2  2.  2.  2  2 


us 


—  ^  Oui  4  Cu2  4  0 U3  4  Cu4  4  Cus  yi  4^V2  4  Cy/44-CyS^. 


where  the  subscripts  1,  2,  3,  4,  5  denote  the  triangle  to  which  the 
notation  refers. 


1 

More  generally,  if  there  are  n  triangles  with  apices  at  0, 


4  (Cuz  4Cv»)hu2  4 . +  (C  un  4  Cy^  )  hUn 

“  (C  Ut^Cyt  +  Cw2  +  Cv2  4  -  -  -  -4  Cyn4Cvn)^o  s  O., 


(3.43 ) 


If  the  common  apex  of  several  triangles  lies  on  an  impermeable 
boundary,  as  shown  in  FIGURE  3.6,  the  finite  difference  equation  is  obtained 
from  equation  (3,43)  by  a  similar  approach  to  that  used  for  obtaining  equation 
(3.42).  An  additional  zone  of  zero  permeability  is  introduced  in  the  triangu¬ 
lar  area  between  the  impermeable  boundaries.  If  these  boundaries  are  located 
next  to  triangles  1  and  n,  equation  (3.43)  reduces  to: 


1  After  this  thesis  was  completed,  a  derivation  of  the  finite  difference  equation 
(3.43)  was  obtained  by  the  writer,  (see  Appendix  F). 


. 
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FIGURE  3.6.  An  apex  of  several  triangles  at  an  impermeable  boundary 


CU|.HU1  -+■  (Cu2+^Vl)^U'2  +  (Cu3  +  - +(^'UO'tC'Vn-i)^un'tCvn.ll 

/  n  *  y  0.44) 

“  (Cut  +  CV|+  CU2  +*  CV£  -b  .  .  -  -b  C  un  +CvnJ  ho  =  O* 

Equations  (3.43)  and  (3.44)  apply  to  cases  where  the  Ou  and  Ov  axes  of  each 
triangle  pass  through  □.  If  some  of  the  triangles  have  their  Ov  and  Ow  axes 
passing  through  0,  the  subscripts  u  and  v  for  those  triangles  are  changed  to 
v  and  w  respectively.  Similarly,  for  any  triangles  having  their  Ow  and  Ou 
axes  passing  through  0,  the  subscripts  u  and  v  are  changed  to  w  and  u. 


3.6  Equations  for  the  residuals 

The  finite  difference  equations  derived  in  SECTIONS  3.3,  3.4  and 
3.5  must  be  satisfied  at  all  appropriate  nodes  for  steady-state  conditions 
to  exist.  It  has  already  been  stated  that  initially  a  guess  is  made  for  the 
head  at  each  node;  hence  the  left-hand  sides  of  the  finite  difference  equations 
are  not  usually  equal  to  zero.  The  following  equations,  which  follow  from 
equations  (3.18),  (3.41),  (3.42),  (3.43)  and  (3.44),  make  allowance  for  this 
difference.  It  appears  in  the  form  of  a  residual,  denoted  by  the  letter  R. 


. 
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For  a  node  inside  a  homogeneous  soil,  the  residual  R  is  given  by 

R  *  Cu.(liu+hu)  +  Cy.(hy  4h") -t-  cw.(hj^  +  hX)  _  ^  (3-45) 

2  (Cu  +  Cv  +  Cw) 

For  a  node  on  a  boundary  between  two  zones  or  triangles  numbered  1  and  2, 


K  = 


£ui-Hu,  +  Cy,.hv,  +•  Cu2-^U2  +  ^V2  +  Cyy|4  Cwe.  {hvV|+  hw2) 


(3.46) 


Cut  +  C  vi  +  Cw.  +  Cue  +  C  V2  +  C 


W2 


where  the  Ow  axes  of  both  zones  lie  along  the  boundary.  For  a  node 
on  an  impermeable  boundary,  but  not  at  an  apex  of  a  triangle, 


Cu-hu  +  Cv- hi  +  C'.i.  h  w  Vi  w 
R  =  _ 2  -  h, 

Cjj  *+■  Cy  +  Cw 


(3.47) 


where  the  Ow  axis  lies  along  the  boundary.  For  a  node  at  the  apex 
of  n  triangles,  numbered  1  to  n ,  not  located  on  a  boundary  of  the  flow  region, 

(C  ut  -*Cyn)-^ui  *  Cyt)hU2.-f - 4  (Cun  ^Vn-|)  ^un  wh^  (3.48) 

Cui  "h^vi  +  Cu2  +  C've.  ■+■ - +  (-un+  ^vn 

where  the  Ou  and  Ov  axes  of  each  triangle  pass  through  the  node. 

For  a  node  at  the  apex  of  n  triangles,  number  1  to  n,  and  also  located  on 
impermeable  boundaries, 

R  =  CU|,  hLH-t  (Cu2-hCvl).hU2-|-(Cu-«>f  Cve)h  — +  (Curt + Cvn-t)  h  un  tCvn^vn  _i 


Cm  +  Cvi  CU2.  +  Cyo  + -  -  +  Cun  +  C 


‘o>(  3.49) 


-vn 


where  the  Ou  and  Ov  axes  of  each  triangle  pass  through  the  node,  and 
the  impermeable  boundaries  are  located  next  to  triangles  1  and  n.  Along  a 
boundary  across  which  water  passes  into  or  out  of  the  flow  region,  the  total 
head  is  constant,  and  can  be  computed  from  the  external  conditions.  Thus  the 
heads  at  all  nodes  located  on  such  a  boundary  are  known  at  the  start  of  the 
analysis,  and,  because  they  are  not  altered  during  the  subsequent  computations, 
the  residual  at  these  nodes  is  zero,  i.e. 
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R  =  O  . 


(3.50) 


The  computation  and  manipulation  of  residuals  at  every  node  in  the 
flow  region,  is  performed  by  the  computer. 


CHAPTER  IV 


COMPUTER  PROGRAM  ASSEMBLY 

4,1  The  functions  of  the  computer  program 

The  computer  proaram  is,  in  effect,  a  list  of  instructions  which 
causes  the  computer  to  perform  certain  operations  and  arrive  at  a  meaning¬ 
ful  result.  The  sequence  of  operations  is  the  same  for  different  sets  of 
data,  until  a  conditional  instruction  is  reached;  the  execution  of  a  con¬ 
ditional  instruction  is  dependent  upon  some  result  obtained  earlier  in  the 
computation.  Proper  use  of  conditional  instructions,  therefore,  enables  the 
same  proqram  to  be  used  for  sets  of  data  requiring  difference  sequences  of 
operations^  . 

for  most  problems  that  call  for  the  use  of  a  computer,  it  is  desir¬ 
able  that  the  data  be  presented  in  its  basic  form,  and  that  the  results  appear 
in  a  readily  usahle  form;  in  this  way,  the  computer  is  doino  as  much  work  as 
possible.  The  basic  data  for  a  seepage  problem  are  the  minimum  and  maximum 
permeabilities  for  each  zone  and  their  directions,  boundary  locations  and 
other  such  items.  The  proqram  is  desianed  so  that  all  the  remaininn  work,  in- 
cludinn  the  division  of  the  zones  into  trianoles,  the  calculation  and  storaae 
of  the  coefficients  'C'  in  enuations  (3.45)  to  (3,49),  and  the  numerical 
approximations,  is  performed  by  the  computer. 

^No  attempt  is  made  herein  to  describe  the  basic  techniques  of  computer  pro- 
aramminq  or  the  functional  components  of  the  computer,  for  this  information, 
reference  should  be  made  to  the  IBM  manual  on  enqineerinq  analysis  for  com¬ 
puters  (  1963  ) , 
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4.2  Description  of  the  computer  program 

A  flow  chart  showinp  dianrammatical 1 y  the  senuence  of  operations 
performed  hv  the  computer  is  presented  in  FIGURE  4.1. 

The  data  is  arranoed  so  that  all  of  the  information  on  one  zone 
is  received  and  the  coefficients  'C*  calculated  for  each  triangle  within  the 
zone  before  going  on  to  the  next  zone ,  (steps  3  to  9 ,  FIGURE  4.1).  In  this 
way,  the  values  of  the  minimum  and  maximum  permeabilities  and  the  direction 
of  the  minimum  permeability  may  be  dispensed  with  as  soon  as  the  coefficients 
' C  *  have  been  calculated.  The  coordinates  of  the  corner  points  of  the  zone 
are  stored  as  apices  of  trianglps,  and  the  same  nronram  statements  are  then 
used  when  considerino  the  next  zone.  The  procedure  specified  for  steps  6, 

7  and  8,  dividino  the  zones  into  trianales,  is  fully  described  in  appendix 
C. 

In  the  data  for  each  zone  the  maximum  desirable  distance  between 
nodes  is  specified.  This  distance  may  vary  from  node  to  node  within  the 
flow  region.  Since  the  sides  of  every  triangle  are  sub-divided  into  the  same 
number  of  segments,  this  number  is  determined  by  the  zone  which  has  the  maxi¬ 
mum  value  of  the  ratio: 

longest  side  of  triangle  in  zone _ 

maximum  distance  between  nodes  for  zone 

Thus  *a  trianole  whose  longest  side  is  10  ft.  and  in  which  the  node  spacino 
is  a  maximum  of  2  ft.  is  the  criterion  rathpr  than  a  60  ft.  triannle  in  which 
a  node  spacing  of  15  ft.  is  satisfactory.  The  decision  to  have  large  node 
soacinas  in  some  zones  and  small  node  spacinns  elsewhere  should  be  made 
cautiously,  because  allowable  errors  arisina  from  thp  coarse  mesh  are  likely 
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(FIGURE  4.1  continued  on  next  pane) 


■ 
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F  rom  10 


(FIGURE  4.1  continued  on  npxt  peqe) 
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From  39 


FIGURE  4.1 


Flow  chart 
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FIGURE  4.2 


Numbering  convention  for  nodes 


to  Rncroach  into  the  fine  mech,  where  they  may  constitute  serious  errors; 
this  effect  is  discussed  in  SECTION  6.1.  because  th*5  number  of  memory  stor- 
aae  locations  in  the  computer  is  limited,  the  program  includes  a  provision  that 
the  number  of  nodes  on  the  sides  of  any  triangle  shall  not  exceed  eleven,  (step 
1 1 ,  FIGURE  4.1) . 

The  arrays  IUU(J2),  1W(J3),  and  I WW ( J 1  )  are  formed  in  step  13  in 
order  that  the  determination  of  the  value  of  the  residual  at  a  node  on  a  tri¬ 
angle  boundary  is  facilitated.  Each  triangle  boundary  lying  on  a  constant  poten¬ 
tial  face  of  the  flow  region  is  identified  by  the  appearance  of  the  number  of 
its  triangle  in  these  arrays.  The  arrays  IUI(J2),  IVI(J3),  and  IWI(Jl)  in 
step  17  perform  a  similar  function  for  the  impermeable  boundaries  of  the  flow 
region . 


The  boundary  conditions  are  completely  defined  and  the  meshes  in 
each  triangle  are  fixed  before  step  19  is  reached.  Any  node  in  the  flow 
region  is  conveniently  identified  by  three  integers;  the  first  integer  is  the 
number  of  the  triangle  in  which  the  node  lies.  Inspection  of  FIGURE  4.2  shows 
that  a  node  lies  on  a  triangle  boundary  if  any  of  the  following  three  conditions 
apply : 


( a 

)  the 

second 

of  the 

three 

integers  is  1 ,  indicating  that 

the  Ov  axis 

the 

boundary  are 

colinear ; 

(b 

)  the 

third  i 

n teger 

is  1  , 

indicating  that  the  Ou  axis  and 

the  boundary 

col 

inear; 

(c 

)  the 

sum  of 

the  second  and  third  integers  is  N+1 ,  where 

N  is  the 

number  of  nodes  on  a  triangle  side,  indicating  that  the  Ow  axis  and  the  boun¬ 


dary  are  colinear 
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Steps  19  to  47  comprise  one  iteration,  in  which  the  residuals  R 
are  computed  for  pverv  node  within  the  flow  renion.  For  any  node  situated  in¬ 
side  a  triangle,  R  is  calculated  usinn  equation  (3.45).  Each  triangle  side 
is  common  to  a  constant  potential  face,  an  impermeable  boundary,  or  the  side 
of  another  trianalp;  thus,  for  a  node  on  a  trianole  side,  the  arrays  IUU(J2), 
etc.,  IUI(J2),  etc.  and  those  .listing  the  coordinates  of  the  aoicps  of  tri¬ 
angles  are  inspected  (steps  36,  37  and  3B )  so  that  the  appropriate  eouation 
for  calculation  R  can  be  selected.  The  procedure  used  for  finding  the  sides 
of  triangles  which  lie  alone  a  common  boundary  is  described  in  appendix  C. 

The  node  at  an  anex  of  a  trianqle  is  oenerally  common  to  other 
trianales,  (FIGURE  3.5).  Each  of  these  triangles  is  searched  for  before 
the  residual  at  the  node  can  bp  computed.  The  node  is  identified  at  first  E)y 
a  narticular  trianqlp,  which  may  be,  for  example,  in  FIGURE  3.5,  trianole  5. 
The  expressions  (Cu5„  h^  +  C^.  h^)  and  (C^  +  Cvc^ )  are  calculated  and 
stored  for  use  later  as  components  of  the  eouation  for  the  residual.  The  pro¬ 
cedure  described  in  an  earlier  paragraph  for  finding  an  ad/joinino  triangle 
is  used  to  determine  the  trianqle  which  is  situated  next  to  trianole  5  nro- 
qressina  in  a  clockwise  sense  about  the  anex,  (step  26).  In  FIGURE  3.5,  tri¬ 
angle  4  is  revealed;  the  expressions  (C  0  h'  +  C  . .  h'  )  and  (C  +  C  ) 
are  calculated  and  added  to  the  expressions  already  obtained.  Eventually, 
every  trianqle  is  revealed  and  thc  rpsidual  is  computed  by  determining  the 
ratio  of  the  sums  of  expressions  obtained  for  each  triangle,  and  subtracting 
from  it  the  existing  value  of  the  head  at  the  node.  If,  while  progressing 
from  triangle  to  trianole  round  the  anex,  a  constant  potential  face  is  encoun¬ 
tered,  the  residual  ic  made  enual  to  zero.  If  an  impermeable  boundary  is  en¬ 
countered,  further  progression  in  a  clockwise  sense  is  prevented;  sn, 
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starting  aaain  at  the  oriqinal  triangle,  any  remainina  triangles  are  re¬ 
vealed  in  counter-clockwise  order,  (steps  20  to  35).  Thus,  in  FIGURE  3.6, 
if  the  node  is  identified  initially  by  trianale  2,  the  next  trianqle  revealed 
is  trianale  1  and  then  the  impermeable  boundary  Ou^  is  encountered;  proceeding 
next  in  a  counter-clockwise  sense  from  trianale  2,  trianale  3  is  revealed. 

After  the  residual  at  a  node  is  calculated,  it  is  multiplied  by  a 
number  qreater  than  one,  called  an  over- relaxation  factor,  and  the  product 
is  added  to  the  current  value  of  the  head  to  obtain  a  new  value,  (step  45). 

The  presence  of  the  over-relaxation  factor  is  discussed  in  SECTION  6.2;  its 
effect  aenerally  is  to  nuicken  the  liouidation  of  the  residuals,  and  so  reduce 
the  number  of  iterations.  If  the  residuals  for  the  nodes  in  an  area  of 
appreciable  size  are  less  than  the  pre-stated  tolerable  value,  but  are  all 
positive  or  all  nerative,  a  larne  error  mav  exist  in  the  calculated  heads  at 
nodes  in  the  centre  of  the  area.  The  use  of  an  over-relaxation  factor  is 
believed  to  overcome  this  effect,  bv  causinn  the  converoence  of  the  head  at 
any  node  to  be  an  oscillatinn  Process.  To  prevent  the  oscillations  becoming 
diveraent,  provision  is  made  for  reducinn  the  over-relaxation  factor  as  soon 
as  ary  increase  occurs  in  the  absolute  value  of  the  residuals  from  one  iteration 
to  the  next,  (step  47). 

Several  sets  of  data  may  be  analysed  in  the  same  run  of  a  nronram. 
After  the  output  for  anv  set  of  data  is  completed,  a  statement  is  read  which 
returns  the  control  to  the  start  of  the  nronram,  and  the  next  set  of  data  is 
input  and  analysed.  The  computer  stops  when  a  dummy  innut  is  received;  this 
innut  consists  of  a  negative  number  or  zero  reoresentino  NZ,  the  number  of 


zones  in  a  flow  renion 


49 


4.3  The  format  of  the  program  and  the  data 

The  oroqram  is  written  in  FORTRAN  language,  which  has  a  format 
close  to  that  of  mathematical  equations  and  uses  certain  words  of  the 
English  lanauaqe  to  describe  each  operation.  The  complete  proqram  is 
presented  in  appendix  D.  Each  statement,  or  line,  is  punched  on  a  stan¬ 
dard  punch  card. 

Built  into  the  proqram  are  a  number  of  statements  which  cause  the 
computer  to  stop  if  the  number  of  items  of  data  is  too  large  or  if  there  is 
a  detectable  error  in  the  data.  Such  an  error  mioht  be  that  the  coordinates 
locating  constant  potential  faces  and  impermeable  boundaries  do  not  agree 
with  the  zone  boundaries.  When  the  procedure  for  selecting  the  triangle 
which  adjoins  a  given  trianale  is  beina  used,  a  set  of  operations  is  repeated 
using  a  subscript  which  is  being  incremented.  Although  this  procedure  is 
desioned  so  that  the  trianqle  is  always  revealed  before  the  subscript  exceeds 
the  total  number  of  triangles  in  the  flow  reqion,  a  statement  preventing  the 
subscript  increasing  indefinitely  is  inserted  into  the  program  as  a  safe¬ 
guard  against  the  computer  wasting  time;  the  statement  causes  the  computer 
to  stop.  Before  any  such  ’unintended  stop’  statement  is  read,  a  messaae 
indicating  the  nature  of  the  trouble  is  passed  through  output. 

In  order  that  each  set  of  data  be  interpreted  meaningfully,  it  is 
arranged  in  a  definite  format.  The  data  is  punched  on  a  series  of  cards  as 
indicated  in  appendix  D.  The  sets  of  data  are  placed  in  order  behind  the 


program  statements. 


CHAPTER  V 


PRESENTATION  OF  WORKED  EXAMPLES 

5,1  Discussion  of  FIGURES  5.1  to  5.6 

Several  problems  for  which  the  computer  program  was  used  are 
shown  in  FIGURES  5.1  to  5.6.  The  largest  residual  remaining  after  the 
final  iteration  is  shown  for  each  case;  this  is  equal  to  or  greater  than 
the  tolerable  maximum,  depending  on  whether  the  specified  number  of 
iterations  was  sufficient  to  allow  the  residuals  to  decrease  to  the  amount 
desired.  The  largest  residual  may  be  regarded  as  an  indicator  of  the 
accuracy  of  the  solution.  In  some  cases,  the  solution  by  other  methods 
is  superimposed  on  the  computer  solution  so  that  a  direct  comparison  may 
be  made. 


In  the  example  of  FIGURE  5.1,  the  flow  region  consists  of  a  homo¬ 
geneous  isotropic  soil  underlying  an  impermeable  structure.  The  flow  is 
confined  by  impermeable  boundaries  at  the  sides  and  the  base,  but  the  dis¬ 
tance  from  the  structure  to  the  impermeable  boundaries  is  large  in  relation 
to  the  width  of  the  structure;  thus  the  flow  pattern  at  all  points  except 
near  the  boundaries  is  similar  to  that  occurring  in  a  semi-infinite  soil, 
for  which  a  mathematically  exact  solution  is  available.  The  complex 
variable  theory  can  be  used  to  show  that  the  equipotential  lines  for  the 
semi-infinite  case  are  a  family  of  hyperbolas  having  common  foci  at  the 
corners  of  the  structure.  Inspection  of  FIGURE  5.1  shows  that. there  is 
close  agreement  between  the  equipotential  lines  obtained  from  the  computer 
solution  and  those  obtained  from  the  complex  variable  theory.  Apart  from  an 


50 


51 

expected  divergence  of  these  lines  near  the  impermeable  boundary,  the 
greatest  difference  between  the  computer  and  complex  variable  solutions  is 
about  one  percent  in  terms  of  potential  head. 

Before  the  data  for  the  example  of  FIGURE  5.1  was  used  by  the 
computer,  the  flow  region  was  arbitarily  divided  into  four  "zones’'  so  that 
the  nodes  were  more  closely  spaced  near  the  structure  than  elsewhere. 

This  procedure  illustrates  how  the  node  spacing  may  be  varied  within  a 
zone  by  simple  menipulation  of  the  data. 

FIGURE  5.2  shows  the  solutions  from  three  different  methods  to 
the  problem  of  flow  round  a  cut-off  curtain  extending  vertically  from  the 
centre  of  an  impermeable  structure.  A  separate  test  has  shown  that  the 
equipotential  lines  as  shown  for  the  electric  analogue  solution  are  almost 
identical  to  those  occurring  in  the  case  of  an  infinitely  long  layer  of  soil. 
Close  agreement  exists  between  the  solutions  obtained  by  the  computer  method 
and  the  conventional  iteration  method  using  a  square  mesh;  the  equipotential 
lines  in  each  case  agree  to  within  about  3  percent  in  terms  of  the  potential 
head  values,  with  the  greatest  difference  occurring  near  the  cut-off  curtain. 
Vertically  below  the  cut-off  curtain  the  solution  using  the  square  mesh  is 
exact,  because  the  symmetry  of  the  flow  region  was  used  in  obtaining  this 
solution.  Although  symmetry  could  have  been  used  in  the  computer  solution, 
it  was  not;  the  error  in  the  solution  near  the  cut-off  curtain  occurs  as  a  re¬ 
sult  of  residuals  not  being  completely  liquidated.  The  electric  analogue 
solution  is  within  10  percent  of  the  iteration  method  using  a  square  mesh, 
with  the  greatest  difference  occurring  near  the  cut-off  curtain.  Experimental 
errors  are  believed  to  account  for  the  difference  between  the  analogue  and  num¬ 


erical  solutions. 
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The  equipoten tial  lines  obtained  by  the  computer  method  for  flow 
through  a  layered  soil  are  shown  in  FIGURE  5.3(a).  The  flow  region  in 
this  example  consists  of  two  anisotropic  zones  and  one  isotropic  zone,  and 
is  confined  at  the  sides  and  the  base  by  impermeable  boundaries.  FIGURE 
5.3(b)  shows  the  equipotential  lines  for  the  same  flow  region  as  FIGURE 
5.3(a),  but  the  data  for  FIGURE  5.3(b)  was  adjusted  so  that  each  zone 
was  arbitarily  divided  into  smaller  "zones".  There  is  a  slight  distinction 
between  the  two  solutions,  attributable  to  the  different  residuals  and  node 
spacings  in  each  case. 

In  FIGURE  5.4,  the  soil  is  homogeneous  and  anisotropic  with  the 
stratification  occurring  at  an  angle  of  25  degrees  to  the  horizontal.  A 
check  was  made  to  verify  the  position  of  the  equipotential  line  representing 
50  units  of  head  by  sketching  the  equivalent  transformed  section  and  assuming 
that  this  line  would  then  be  approximately  perpendicular  to  the  base  of  the 
structure.  Using  the  transformed  section,  the  value  obtained  for  the  angle 
between  this  line  and  the  horizontal  in  the  flow  region  before  transformation 
was  72 degrees;  the  value  of  the  corresponding  angle  in  FIGURE  5.4  is  71y 
degrees • 

FIGURES  5.5  and  5.6  are  illustrations  of  the  applicability  of  the 
computer  program  to  flow  regions  consisting  of  irregularly  shaped  zones.  A 
partial  check  of  the  solutions  was  made  by  verifying  the  continuity  of  flow 
condition  normal  to  a  zone  boundary.  The  normal  velocity  components  at 
eight  randomly  chosen  points  on  boundaries  were  computed  using  the  heads 
shown;  the  computed  values  on  both  sides  of  a  boundary  were  then  compared 
for  equality.  A  difference  of  5  to  8  percent  was  evident  in  all  cases,  and 
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was  probably  caused  by  the  fact  that  the  velocity  component  was  an  average 
value  for  some  length  normal  to  the  boundary.  These  differences  were 
neglected  for  the  purpose  of  verifying  the  continuity  condition. 
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NUMBER  OF  ITERATIONS:  150 
LARGEST  RESIDUAL:  0025 
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NUMBER  OF  ITERATIONS:  120 
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ISOTROPIC 

LARGEST  RESIDUAL.  0  013 

EXCEPT  THAT  EACH  ZONE  IS  DIVIDED  INTO 
SMALLER  "ZONES"  ALONG  DOTTED  LINES 


FIGURE  5-3  FLOW  UNDER  A 


STRUCTURE 
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FIGURE  5-4  FLOW  THROUGH  A  HOMOGENEOUS  ANISOTROPIC  SOIL 

UNDERLAIN  BY  AN  IMPERMEABLE  STRATUM 
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CHAPTER  VI 


DI5CUS5I0N 


6,1  The  approximation  in  the  finite  difference  expressions 

The  finite  difference  equations  derived  in  CHAPTER  III  are 


approximations  to  the  true  continuity  equation  for  flow  through  an 
elemental  area  in  a  plane.  They  are  approximate  because  the  expression 
for  the  second  differential  of  head  with  respect  to  distance  along  an 
axis  in  terms  of  head  values  at  nodes  on  that  axis  is  not  exact.  For 


example,  along  the  Ou  axis  this  expression  is  stated  as  equation  (3,14a). 


(3.14a) 


The  terms  in  the  Taylor*s  series  expansions  for  h^  and  h^  which 
are  omitted  when  setting  up  equation  (3.14a)  are 


uuuu 


+  terms  of  the  order 


etc . 


The  error  in  the  finite  difference  equation  for  the  head  at  a  node  not  lying 


on  a  boundary  is  therefore  made  up  of  three  parts,  each  proportional  to  the 
square  of  the  node  spacing  along  the  three  axes.  Ignoring  temporarily  the 
effect- of  the  fourth  derivatives  of  h  in  any  direction,  the  largest  node 
spacing  along  any  axis  through  a  point  controls  the  accuracy  with  which  the 
head  at  that  point  can  be  computed.  After  several  iterations,  an  error  in 
the  head  at  one  node  is  manifested  in  the  head  at  nodes  inside  an  area  of 
appreciable  size.  Therefore,  not  counting  the  difference  in  permeability 
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from  one  to  another,  or  differences  in  the  fourth  derivative  of  head 

throughout  the  flow  region,  the  accuracy  of  the  head  at  all  points  is  in 
general  controlled  by  the  largest  node  spacing  throughout  the  region,, 

This  factor  should  be  considered  when  the  choice  of  node  spacings  is  being 

made,  but  it  need  not  be  a  deterrent  to  the  choice  of  large  and  small  node 

» 

spacings  in  a  flow  region,  if  it  is  remembered  that  the  accuracy  in  the 
areas  where  the  nodes  are  closely  spaced  may  be  reduced  by  the  existence 
of  the  large  node  spacings  elsewhere.  Large  and  small  node  spacinas  are 
chosen  for  the  example  shown  in  flLURE  5.1,  and  thus  the  accuracy  of  the 
heads  at  points  near  the  structure  are  probably  reduced  by  the  large  node 
spacings  near  the  impermeable  boundaries;  however,  any  loss  of  accuracy 
appears  to  ba  slight  in  terms  of  practical  values,  and  is  therefore  out¬ 
weighed  by  the  fact  that  closely  spaced  nodes  help  in  obtaining  a  better 
idea  of  the  flow  pattern  near  the  structure. 

The  error  in  the  finite  difference  equation  is  proportional  to  th 
fourth  derivative  of  h.  Throughout  most  of  a  flow  region,  this  value  is 
negligible,  but  there  might,  in  some  cases,  be  points  where  this  value 
approaches  infinity.  5uch  points  exist  where  the  superficial  velocity  is 
theoretically  infinite.  Thus  at  these  points  finite  difference  expressions 
such  as  equation  (3.14a)  are  incorrect;  however,  a  small  distance  away  the 
expressions  are  correct,  and  for  this  reason  their  use  at  every  point  in 

i 

a  flow  region  appears  justified. 

6.2  Over-relaxation 

Because  the  methods  of  numerical  approximation  were  devised  with 
the  intention  that  solutions  be  obtained  by  hand,  it  is  frequently  recom¬ 
mended  that  the  estimates  for  the  heads  initially  assigned  to  each  node 
should  be  as  accurate  as  possible,  so  that  the  amount  of  subsequent  work  is 
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minimised.  As  a  result,  the  negative  and  positive  residuals  are  randomly 
distributed  throughout  the  flow  region.  However,  when  the  solution  is  being 
obtained  by  computer,  it  is  more  convenient  to  assign  the  same  initial 
value  of  head  to  every  node  except  those  lying  on  a  constant  potential  face, 
even  though  this  value  may  be  an  obviously  poor  guess.  The  number  of 
iterations  required  to  obtain  a  solution  is  then  generally  greater  than  if 
a  good  guess  were  made  and  in  order  to  reduce  this  number,  over-relaxation 
is  used.  The  residual  at  any  node  is  multiplied  by  an  over-relaxation  factor 
before  adding  it  to  the  current  value  of  the  head. 

The  following  is  a  hypothetical  case  to  illustrate  the  effect  of 
over-relaxation .  Consider  a  node  at  which  the  true  value  of  the  head  is 
75,  and  the  value  initially  assigned  is  50.  The  values  of  the  head  after 
successive  iterations  in  which  over-relaxation  is  not  used  are,  in  order, 

50,  65,  70,  73,  74,  74.5,  74.8,  74.9,  74.95  . 

The  values  after  iterations  in  which  over-relaxation  is  used  are 

50,  70,  76,  74.5,  75.1 ,  74.95  . 

It  is  seen  that  the  heads  converge  quickly  with  the  use  of  over-relaxation; 
and  that  the  convergence  is  an  oscillating  process.  It  is  important  that  there 
be  a  random  distribution  of  negative  and  positive  residuals  at  all  heads  in  the 
flow  region.  Since  the  initial  values  of  the  head  at  all  nodes  is  the  same, 
the  residuals  initially  are  all  positive  in  some  areas  and  all  negative  in 
others,  and  it  is  the  process  of  over-relaxation  that  disperses  the  signs 
of  the  residuals  after  sufficient  iterations  have  been  made  for  the  convergence 
to  oscillate.  It  therefore  appears  necessary  to  over-relax  the  residuals  for 
two  reasons,  to  quicken  the  convergence  and  to  ensure  a  random  distribution 
of  negative  and  positive  residuals. 
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The  writer  was  unable  to  find  a  logical  procedure  for  determining 
the  over-relaxation  factor.  At  first  the  value  of  1*1  was  chosen,  but  the 
effect  on  the  rate  of  liquidation  of  the  residuals  for  a  number  of  examples 
was  found  to  be  negligible.  The  values  1 .5  and  1 .7  were  tried,  and  finally 
1.6  was  decided  upon.  For  a  flow  region  similar  to  that  shown  in  FIGURE 
5.1,  the  number  of  iterations  required  to  obtain  a  maximum  residual  of  0.010 
was  over  1000  when  no  over-relaxation  was  used;  when  an  over-relaxation 
factor  of  1 .0  was  used,  the  number  of  iterations  for  the  same  maximum  resid¬ 
ual  was  109.  In  some  examples,  the  oscillations  of  the  residual  at  certain 
nodes  becomes  divergent  with  the  use  of  over-relaxation.  This  situation  is 
believed  to  be  caused  by  the  interaction  of  in-phase  oscillations  at  adjacent 
nodes.  The  writer's  method  for  preventing  the  growth  of  divergent  oscilla¬ 
tions  was  to  reduce  the  over-relaxation  factor  immediately  it  arises.  The 
amount  by  which  the  over-relaxation  factor  exceeds  one  is  multiplied  by  0.8; 
thus  if  the  tendency  to  diverge  occurs  three  times  before  the  final  iteration, 
the  over-relaxation  factor  becomes  equal  to 
1  +  0.6.  (0.8)^t  or  1.41. 

Although  the  procedure  for  determining  the  over-relaxation  factor 
is  somewhat  empirical,  it  appears  to  perform  its  function  satisfactorily. 

6.3  Methods  of  manipulating  the  residuals 

The  three  following  methods  of  manipulating  the  residuals  were 
considered  before  the  f rogram  was  assembled  in  its  present  form:- 

(a)  After  the  head  at  a  node  is  adjusted  by  the  amount  of  the  residual 
multiplied  by  an  over-relaxation  factor,  the  residual  at  that  node  is  not 
required  any  further,  and  when  the  residual  for  the  next  node  is  calculated, 
it  occupies  the  same  storage  space  in  the  computer. 
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(b)  The  residual  is  computed  and  stored  as  an  element  of  a  three- 
dimensional  array.  After  the  residuals  at  all  nodes  are  computed,  the 
heads  are  adjusted. 

(c)  The  residual  is  computed  and  stored  as  an  element  of  a  three- 
dimensional  array.  After  all  the  residuals  are  computed,  they  are  relaxed 
by  adjusting  the  head  at  each  node  and  the  residuals  at  surrounding  nodes. 

In  the  relaxation  procedure  by  hand  using  a  square  mesh,  method 
(c)  is  used  most  frequently  because  much  of  the  arithmetical  work  involves 
only  small  numbers.  It  is  possible  to  see  at  a  glance  the  locations  of  the 
largest  residuals  in  the  flow  region  and  reduce  these  a  number  of  times  be¬ 
fore  considering  other  nodes.  This  method  is  not  adopted  in  the  computer 
solution  because  a  large  number  of  memory  storage  locations  would  be  required 
for  the  residuals,  and  the  determination  of  the  adjustments  to  be  made  to 
the  residuals  at  surrounding  nodes  would  be  a  lengthy  operation.  Method 
(a)  has  the  advantages  that  it  does  not  require  a  large  number  of  memory 
storage  locations  for  the  residuals,  and  the  number  of  iterations  is  less 
than  would  occur  in  method  (b)  because  the  heads  are  adjusted  after  the 
calculation  of  each  residual.  For  these  reasons  method  (a)  is  adopted  in 
the  computer  program. 

6.4  Validity  of  the  finite  difference  equation  :  point  at  the  apex  of  a 

triangle 

In  SECTION  3.5,  a  finite  difference  equation  for  a  point  at  the 
apex  of  a  triangle  (equation  3.43)  was  deduced  but  not  proved.  The  date  for 
the  example  of  FIGURE  5.3(b)  was  designed  so  that  this  equation  could  be 
tested.  FIGURES  5.3(a)  and  (b)  both  represent  the  same  flow  region,  but 
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FIGURE  5.3(b)  is  divided  into  a  greater  number  of  triangles,  and  con¬ 
sequently  equation  (3.43)  i9  used  many  more  times  in  the  analysis  of 
FIGURE  5.3(b)  than  it  is  in  FIGURE  5.3(a).  The  fact  that  the  solutions 
from  both  sets  of  data  agree  closely  is  evidence  for  the  soundness  of  the 
equation • 

6.5  The  accuracy  of  the  solutions 

The  close  agreement  between  the  solutions  from  the  computer  and 
other  methods  for  the  examples  of  FIGURES  5.1  to  5.3  substantiates  the 
validity  of  the  computer  program.  The  accuracy  of  the  computer  solutions 
could  be  increased  by  specifying  a  smaller  tolerable  residual  and  a  greater 
allowable  number  of  iterations.  For  practical  purposes,  accuracy  greater 
than  that  shown  in  FIGURE5  5.1  to  5.6  is  hardly  warranted  and  a  tolerable 
residual  equal  to  about  0.02  percent  of  the  total  head  difference  appears 
satisfactory  for  flow  regions  having  the  complexity  shown. 

Another  factor  to  be  considered  when  deciding  the  allowable  number 
of  iterations  is  the  cost  of  computer  time  in  relation  to  the  accuracy  obtained. 
For  the  flow  region  shown  in  FIGURE  5.1,  the  computation  time  required  by 
the  IBM  7040  computer  used  by  the  writer  was  3  minutes;  for  the  more  complex 
flow  region  of  FIGURE  5.6,  the  time  was  nearly  19  minutes.  Until  more  flow 
regions  have  been  analysed  using  the  program,  it  is  not  passible  to  present  a 
simple  relation  for  estimating  computer  time;  however,  it  is  likely  that  an 
appreciable  proportion  of  this  time  is  spent  in  searching  for  the  sides  of 
triangles  which  lie  along  a  common  boundary. 


CHAPTER  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 

7.1  Conclusions 

In  existing  methods  of  analysing  the  flow  of  water  through  non- 
homogeneous  anisotropic  soils,  simplifying  approximations  in  the  data  are 
necessary  in  order  that  a  solution  may  be  found  within  a  reasonable  time. 
Although  these  approximations,  if  made  with  discretion,  may  be  small  in 
terms  of  the  effect  on  the  final  result,  any  attempt  to  devise  a  method  of 
analysis  in  which  these  approximations  are  not  required  is  justified  for 
engineering  purposes,  if  the  procedure  for  solution  involves  no  cumbersome 
operations  and  can  be  completed  quickly.  The  writer  concludes  that  such  a 
method  has  been  devised  in  this  thesis. 

The  basic  two-dimensional  flow  equation  is  developed  in  finite 
difference  form  in  terms  of  three  axes  in  a  plane.  The  principles  of 
numerical  analysis  are  then  used  to  solve  this  equation  at  nodes  of  tri¬ 
angular  meshes  superimposed  on  the  flow  region.  Each  triangular  mesh  is 
formed  by  arranging  that  certain  axes  lie  along  boundaries  between  zones  of 
soil  having  different  permeability,  as  explained  in  SECTION  2.4.  Conventional 
methods  of  numerical  analysis  for  flow  problems  use  a  square  mesh,  (Shaw  and 
Southwell,  1941),  but  since  the  axes  of  a  square  mesh  cannot  usually  be 
arranged  to  coincide  with  the  boundaries  of  zones  within  a  flow  region,  the 
finite  difference  equations  for  nodes  situated  less  than  one  node  spacing 
from  the  boundaries  are  awkward  to  set  up.  By  adjusting  the  mesh  configuration 
to  suit  the  boundaries,  as  is  done  in  the  triangular  mesh,  the  problem  of 
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such  nodes  does  not  arise  and  the  formation  of  the  finite  difference 
equation  for  any  node  in  the  flow  region  becomes  systematic. 

An  important  advantage  of  the  method  is  that  it  is  programmed  for 
solution  by  computer.  Due  to  the  versatility  of  the  triangular  mesh  con¬ 
figuration,  the  program  can  be  used  for  the  analyses  of  flow  regions  of  any 
desired  complexity,  subject  to  the  limitations  imposed  by  the  size  of  the 
computer,  which  are  listed  in'  SECTION  1.3. 

The  program  is  limited  to  cases  where  no  phreatic  line  is  present, 
and  cannot  therefore  be  used  for  analysing  the  flow  of  water  through  earth 
dams  unless  the  position  of  the  phreatic  line  is  assumed.  However,  the 
writer  believes  that  the  program  could  be  extended  to  cover  cases  which 
include  a  phreatic  line,  and  his  proposals  in  this  respect  are  outlined  in 
SECTION  7.2. 

The  validity  of  the  computer  program  is  substantiated  by  the  close 
agreement  shown  in  the  examples  of  FIGURES  5.1  to  5.3  between  the  solutions 
obtained  from  the  computer  and  from  other  methods.  The  general  applicability 
of  the  program  is  illustrated  by  the  flow  regions  in  FIGURES  5.5  to  5.6  for 
which  a  solution  by  any  other  reasoned  procedure  would  be  impractical. 

7.2  Recommendations 

The  computer  program  as  presented  in  this  thesis  is  by  no  means 
utterly  efficient  in  the  way  the  procedure  for  analysis  is  specified  to  the 
computer.  This  situation  is  largely  due  to  the  fact  that  the  writer  has 
concentrated  on  getting  the  method  working,  rather  than  on  procedural  details 
within  the  program  itself.  For  example,  the  program  could  be  improved  by 
altering  the  steps  involved  in  searching  for  a  triangle  side  which  adjoins 
a  given  line  so  that  after  the  search  has  been  made  during  the  first  iteration, 


JjuuTitf. 
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the  number  identifying  the  triangle  side  is  replaced  in  a  memory  storage 
location  in  such  a  way  that  it  may  be  obtained  without  a  search  in  the 
second  and  subsequent  iterations.  No  doubt  further  improvements  can 
be  made  which  would  reduce  the  length  of  computation  time,  and  it  is 
suggested  that  some  effort  be  made  in  this  direction  before  any  work  is 
done  in  extending  the  program  in  any  way  outlined  below. 

Along  boundaries  across  which  wgter  leaves  the  flow  region  it  is 
important  that  the  hydraulic  gradient  be  less  than  that  which  would  cause 
removal  of  the  soil  from  the  surface,  i.e.,  piping.  In  this  connection,  it 
is  recommended  that  the  computer  method  be  checked  to  determine  if  sufficient 
accuracy  is  obtained  in  the  vicinity  of  these  boundaries.  In  a  strict  sense, 
computation  of  the  hydraulic  gradient  requires  that  the  direction  of  flow 
be  known;  this  direction  is  perpendicular  to  the  equipotential  line  at  any 
point  in  an  isotropic  soil,  but  this  rule  does  not  exist  in  an  anisotropic 
soil.  The  direction  of  flow  can  be  obtained  from  a  knowledge  of  the  flow  lines, 
and  it  is  suggested  that  these  may  be  determined  by  computing  the  velocity 
components  in  the  principal  permeability  directions  at  nodes  situated  where 
the  flow  lines  are  required. 

A  natural  extension  of  this  thesis  would  be  an  adaptation  of  the 
program  to  cover  flow  regions  which  have  a  phreatic  line  as  a  boundary.  The 
writer  suggests  that  this  adaptation  could  be  made  by  initially  estimating 
the  location  of  the  phreatic  line,  in  much  the  same  way  that  an  initial  guess 
was  made  for  the  head  at  each  node,  and  then  treating  that  line  as  if  it  were 
an  impermeable  boundary;  if  the  line  is  curved  it  is  approximated  to  a  number 
of  straight  lines.  After  a  sufficient  number  of  iterations  has  been  made  to 
cause  residuals  at  nodes  to  become  tolerably  small,  the  trial  phreatic  line  is 
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inspected  for  the  condition  that  the  change  in  total  head  between  any  two 
points  along  its  length  iB  equal  to  the  difference  in  vertical  elevation  between 
those  points.  The  location  of  the  phreatic  line  is  altered  until  this  con¬ 
dition  is  satisfied  within  tolerable  limits.  With  each  alteration  the  mesh 
configuration  is  adjusted  to  suit  the  new  alignment  of  the  phreatic  line 
and  new  values  for  the  heads  at  nodes  are  determined.  In  effect,  the  procedure 
for  analysing  flow  regions  which  do  not  include  a  phreatic  line  is  repeated 
until  the  location  of  the  phreatic  line,  temporarily  regarded  as  an  imperm¬ 
eable  boundary,  is  determined  correctly  through  successive  approximations. 

The  presence  of  a  seepage  surface  may  be  regarded  for  the  purposes  of  the 
analysis  as  a  boundary  between  two  zones,  one  of  which  is  infinitely  permeable 
and  equations  (3.46)  to  (3.49)  adjusted  to  suit  this  condition. 

It  is  recommended  that  the  concept  of  a  triangular  mesh  be  used 
as  a  basis  for  investigating  steady-state  three-dimensional  flow  conditions 
in  non-homogeneous  soils.  Each  zone  could  be  approximated  to  a  solid  figure 
whose  sides  consist  of  irregular  triangles;  the  zone  can  then  be  divided  into 
a  number  of  tetrahedra.  The  finite  difference  equation  for  a  node  may  be 
expressed  in  terms  of  the  six  axes  which  are  parallel  to  the  sides  of  each 
tetrahedron.  The  positions  of  nodes  can  then  be  classified  under  four 
headings,  as  follows:  (a)  inside  a  tetrahedron,  (b)  on  a  side  of  a  tetrahedron, 
(c)  on  an  edge  of  a  tetrahedron,  nr  (d)  at  an  apex  of  a  tetrahedron;  the  pro¬ 
cedure  for  solution  is  then  likely  to  be  similar  to  that  used  for  two¬ 


dimensional  flow. 
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APPENDIX  A 


CONTINUITY  EQUATION  :  DERIVATION 

Consider  an  arbitary  point  0  in  an  anisotropic  homogeneous  soil 
in  which  flow  occurs  in  two-dimensions  only',  defined  by  the  plane  of  Ox  and 
Oy.  Let  Ox  and  Oy  be  the  directions  of  principal  permeabilities,  and  vx,y 
the  components  of  (superficial)  velocity  at  0,  (FIGURE  A 1  ) . 

The  velocity  components  at  P  on  the  circumference  of  the  elemental 
circle  centred  on  0,  radius  As,  are  equal  to  the  velocity  components  at  0 
plus  the  rate  of  change  of  the  velocity  components  in  the  s-  direction  multi¬ 
plied  by  As,  i.e.l 

velocity  component  at  P  in  x-  direction  =  vy  +  (vy)c;.As,  and 
velocity  component  at  P  in  y-  direction  =  v  +  (v  )  .  As • 

y  y  s 

Therefore  the  flow  q  across  PP*  is  equal  to  the  radial  velocity  component,  mul¬ 
tiplied  by  the  arc  length  PP'.  Resolving  the  above  velocity  components  in 
the  radial  direction, 

^Partial  derivatives  are  denoted  herein  as  follows,  e.g: 

(h)  ,  which  represents  jlll,  or  the  partial  derivative  of  h  with  respect 

u  du 

to  u,  and 

d?h 

(h)  ,  representing  _ ,  or  the  second  partial  derivative  of  h  with 

uv  du .  dv 

respect  to  u  and  v. 

A2 


FIGURE  A1 . 


Elemental  circle. 


(\  =  {(Vfc  +  •  As)  C©Scrf  4-  (vy  4  (vy)s  .  As)  Sin  c<  }  As  .  Arf 


{v^Cos*  +  Vy  Sin  oL  ]”  As.  G* 

{  (Vx)s  <*  +  (Vy)s  Sl'no<}  .  (As)2.  . 


”  (3.1) 


But,  by  Darcy's  lew,  v  =  -k  (h)  and  v  =  -k  (h)  where  k  .  are  the  permea- 

x  xx  y  yy  X  y 

bilities  in  the  x-  and  y-  directions  respectively.  Therefore,  differentiating 


the  expressions  for  v  and  v  , 

x  y 

n 

V) 

k-X  C^OxS  , 

(3.2a) 

X-N 

11 

-ky  (h)yS  . 

(3.2b) 

Substituting  equations  (3.2a)  and  (3.2b)  in  equation  (3.1), 


%  -  [Voc  COS  Crf  4  Vy  Sin  erf  j  .  As  *  <5* 

{  k* (h)xs  005 *  4  ^  (h)ys  Sin  *}  .  (As)2  &rf  . 


(3.3) 
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It  is  now  necessary  to  express  the  mixed  second  partial  derivatives 

(h)  .  (h)  in  terms  of  the  second  partial  derivatives  of  h  with  respect  to 

xs  ys 

distances  measured  along  the  Ou,  Ov  and  Ow  axes.  (i.e. , (h).  (h)  .  .  (h)  ), 

^  '  *  uu  vv  ww 

where  these  axes  are  parallel  to  the  sides  of  the  triangular  mesh.  Two  axes 
(albeit  skew)  are  sufficient  to  define  a  point  in  a  plane;  the  Ou  and  Ow  axes 
are  arbitarily  chosen  for  this  purpose.  For  a  change  from  the  point  0  to  the 
point  P  (by  As),  the  u-  dimension  changes  by  Au>  where 


As.  Sin  (ot  w  -ot)  . 
Sin  {ot  vj  -tku) 


this  equation  is  obtained  by  use  of  the  sine  rule  on  triangle  OPW 


•  i  i 


(FIGURE  A2 )  • 


In  the  limit,  Au  becomes  ^u  and  As  becomes  As,  so  that 


Sin  y/  -of*) 

Sin  (diyj-oi  u)  5 


(3.4a) 


Similarly,  it  can  be  shown  that 


(w)5  -  -  Sin  (otu-o<) 

Sin  (otw-otu) 


(3.4b) 


and 


(U)*.  =  Sl”  **  (3.4c),  (». 

Sin  (otw-otu)  9 


(3.4d), 

s,n  0**1  -°<u)  * 


(u)r  =  -  006 


Sin  ? 


(  3 . 4e  ) ,  (W)u  = 


_  ( 3 . 4f ) , 

5in  (*w-o*u)  ’ 


(u)v  *  S'n  Ow 

Sin  (<*  w  -  oij  7 


(w)^  =  -  1  n 


(3.4h) 


(3 ,4q ) , 


S  i  n  v/  ~  o*  u  ) 


FIGURE  A 2. 


Continuing  to  work  in  terms  of  the  Ou  and  Ow  axes,  any  distance 
along  Ox,  Oy,  Os  or  Ov  may  be  regarded  as  a  function  of  u  and  w  only;  therefore 

(h)s  =  (h)u-(u)s  +  (h)w.(w)s.  (3.5) 

Substituting  equations  (3.4a)  and  (3.4b)  in  equation  (3.5), 


(h)s  =  (h)u..an  (<„•■')  _  (h)w  .  Sm  . 

Si  n  (d  ~cL m.) 

Differentiating  both  sides  of  equation  (3.6)  with  respect  to  x, 


(3.6) 


Sip  (o*W  -oQ 

Sin  (o<w  -otu) 


X 


Sir\ 
5  m 


(*u  ~<*)  j 

(cXw  -c<u)/x 


9 
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(h)*s  =  f(h),A,  .('uW..Slri  few-*)  4  ((h)u')w  •  (w)x. 


5m  (dw-oiu) 


Sin  (cL  -Ol^ 


>.(3.7) 


_  ((h^  .  (A.  s;«(«u-«) 


Sin  (olw-<A) 


iin  (o<y yj  -O^u) 


Substituting  equations  (3.4c)  and  (3.4d)  in  equation  (3.7),  and  rearranging, 


(h)*5  ®  (h)uu  .  $'no<W>  Sm  +  (la)v/w  .  Sin  c<u.  51  n  (otu 

Si  na  (« w  -  o<u)  Si  n2  (c*w  -ctf  u) 


>(  3 . 8  a  ) 


-(h)uw.  3m  sin  (p<w -op  +  sin  t*w  .  sin  (o(u  -  op 

Sin2  («w-o<u) 

Similarly,  using  equations  (3.4e)  and  (3.4f),  it  can  be  shown  that 

(h)ys  ~  ~  -  CoS(^W  .  Sin  (ols^-od  -  (H)wW  •  ofu.Smfeu-*) 

si *n2  (t**  - <xu)  Sin1  w  -*u) 

.  fu\  Cos  <*u.  sun  (&<w-oO  +  c*w  .  sin  ( oLu-ot. )  . 

Tlr  /uw‘  - - ;— 5-7 - : - 

Sm  (o<w  -  o^uJ 

Using  equations  (3.4g)  and  (3.4h),  it  can  be  shown  that 


>(3.0b) 


(h)w  =  (h)uu.  6ma(»<»v  -«A  4  (h)ww.  Slf,a^u 

<$\r\'2-  (o(w -ofu)  (j?(w -t/u) 

-2(h)uw-  6|n  (^w-c/y).  Sm  (p/u-'cx'v)  # 

Sin2  -c*u) 


>  (3.8c) 


Equation  (3.8c)  may  be  used  to  eliminate  (h)  from  equations 


(3.8a)  and  (3.8b),  resulting  in 
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(h)»s  =  (  .  Sln  °»v  sin  Sin  «w.sin  (ot-oty) 

2  sin  s,n  (oC  u  wt^w) 


+  (h)vv.  5'n  <*v/«  Sin  (ot  -<Xg)  4-  s'un  c<u .  sun  fa-o<w)  (3.9a) 

2  sin  (o^-tfy).  Si n(o(v-<yu) 

-4-(h)ww.  S<n  <*u,  -  sin  (of  -<*v)  4*  sir»  <*v  .S>o  (ot-Uu)  > 

2  Sun  (o/u  —  ol  w)*  S’f'' 


and 

(h)  =  -  (h)uu .  005  *v>  s<n  (°<-o<w)  -+-  Coe  cxw.  sinfc-^y) 


2  Sin  (o(v-<Xu).  3in(ofu-0t w) 


-(h)vw.  ^  S,n  6<  -0<u)  +  <*u  -  S'ynfa-o/vj) 

o  -  •  .  .  \  «•  /•  .  .  \ 


2.Stn  (c<w-{^y).  SirN  (o^v-^u) 

-  (h)^,  C^Q/U>  sin  Gx-qQ  *  Ccs^y  ,  sin(o(^u) 

2  Sin  Sfn(o^vv/~^v)  > 

Equations  (3.9a)  and  (3.9b)  may  be  written  as  follows: 

(h)xs  =  y1  (h)MM  .  SiV>6*-°<p)  +  sin  Up,  Sinfr-tf^)  ? 

. .  ^ — 1  2  sin  (c*  wi  -o^m!  sin 


►  (3.9b) 


(3.9c) 


M=u,VjW 
and 


M  =  y  (h)MM.  coa  dN.  sin  ^-«P)  4-  toi«p  .  Sin(ot-otN)  (3  9d) 

WyS  ^  2  («N-olM).sin  (oiM-o((>) 


M?u;v,w 


where  M,  N  and  P  are,  respectively,  u,  v,  w, 

and  v,  w,  u, 
and  w,  u,  v. 


Substituting  equations  (3.9c)  and  (3.9d)  in  equation  (3.3),  the 


flow  q  across  PP*  is  given  by 


A0 


=  {*Vx  cos  o<  +  Vy  sin  o<|.  A s,  5o< 

—  J  k  ^ ,  COS  o<  Y  ^l^^N>Sin(^-Q<p)'t~S<n^F»5|n(0<-o<M)  |  & 


1  s'in  (o<N-o(M).  sin  (d^-oL?) 


4, 


« 

ky.  sin«V (h)MM.  Ccso(N.s;n  (u-olp)  +  Costir.S in(<x-Vig)  ,  (^2.Sv 

2  sin  sin  (o(M-<Kp)  / 


*  { 


vx  COS  ot  +  Vy  s\n  As. 


-JV(k)  f(kx4-k^)sio(o(N^p)ot^Si^-2k^Sin^N^^o(Fcaso<~2kyGos«<K<»»MfSin,^| 

^  2  Sin  sin  (o^M-o(p)  ) 


<3.10) 


Continuity  requires  that  the  net  outflow  from  the  circle  is  zero. 
This  condition  may  be  stated  algebraicly  as 

2tt 

q  dd  =  O  . 


f 


If  0  is  a  point  inside  a  homogeneous  soil,  then  k^  and  are  constant  values, 

Thus,  integrating  equation  (3.10)  between  0  and  2tt  , 

/>  2tT  r^TT  ..  ^ 

c^doi  =  I  jvxGo->o<  4  Vy  Sind  l.  .  doC 

-  f  r f^xtkg).sir\(<<w40o>6t( S>not-2  S>n <*>&*<* -2kUCoSo{^toSo<pSi nfel  / Ac? 

Jo  l  MML  2  Sin  (c*n-oCm.).  (o<*-c<p) 


-  TT 


(A*)4  ^  (h) 


MM 


k  x  Sin  sin  oiF  4*  ky  c^s  Cos o(f> 

Sin  sin 


Therefore,  for  continuity, 


E/|  \  kx  Slrt  o<N  •  d«"  <*P  +  ky  Ce«oZN.0=6oip 

l  n J  .  - — - t - r - r-*—. - - -  -  O.  (3.11) 

V  Sin  (<*M  -<**,).  Sin  (<** -o<f). 


M*u,v,w 


Equation  (3.11)  is  the  continuity  equation  of  flow,  in  terms  of 


three  axes  in  a  plane. 
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FINITE  DIFFERENCE  EQUATIONS  s  DERIVATION 


Point  inside  a  homogeneous  soil 

Consider  a  point  0  at  a  node  in  a  triangular  mesh  superimposed 
on  a  homogeneous  anisotropic  soil.  The  Ou,  Ov  and  Ow  axes  are  at  angles 
c(u  ,  dLy  and  ^measured  in  a  counter-clockwise  sense  from  the  direction  Ox 
of  minimum  permeability.  The  angle  between  Ox  and  some  datum  direction 
(arbitarily  taken  as  horizontal)  is  «a .  The  points  U',  V'  and  W'  are  one 
mesh  length  from  0  in  a  positive  direction  along  the  Ou,  Ov  and  Ow  axes 
respectively;  U" ,  V"  and  W"  are  one  mesh  length  in  a  negative  direction 
along  the  same  axes.  Let  A-)U,  A^v,  w  be  the  mesh  lengths  along  the  Ou, 

Ov  and  Ow  axes  respectively,  (FIGURE  Bl). 

Denoting  the  head  at  0  by  hQ,  at  U'  by  h^,  and  at  U"  by  h”f  the 
variation  of  head  along  Ou  may  be  expressed  by  Taylor's  theorem  as  follows: 


hi*  =  h0  +  A,u.(h)u  +  (4ll^.(h)uu  +  (A!d).(H)i 

21  3l 


uuu 


%i‘oo 


UUUU 


(3.12a) 


and 


K  -  K  -  A,u.(h)u+(4^\(h)uu_(^’(h)uuu  +  (A^VUuuua+...^  (3. :2b) 

c. .  3  •  • 


Adding  equations  (3.12a)  and  (3.12b)  gives 

(h)uu  = 


(A,u) 


uuu 


(3.13) 
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y 


Lt  axis 
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FIGURE  B1 


By  choosing  a  mesh  size  which  is  small,  terms  containing  A^j  u  raised  to  the 
power  of  2  or  greater  may  be  neglected,  and  the  second  derivative  of  head 
with  respect  to  distance  along  Ou  is  given  approximately  by 


(h)uu  =  hk  *  h»  ~2h°  . 

(A,u)2 

Similarly,  with  respect  to  distances  along  Ov  and  Ow, 

(h)w  =  h;  4-  k  -  . 

(A,  v)2 

and 


hw  +  hw  Zh0 

(k ,  w)2 


(3.14c) 


Now,  referring  to  FIGURE  B1  ,  the  length  U’V*  =  length  OW  =  Aiw. 


Therefore,  using  the  sine  rule  on  triangle  OU'V', 


Bd 


A  |  U 

S\r\  (  yj  —  °^v) 


A.  v 


||  w 


5m  (c<  w  “  of  u) 


Sin  (ofw  -  ot  u) 


Denoting  each  of  the  above  equalities  by  F ,  and  substituting  for  A^u,  A^v, 

A-|  w  in  equations  (3.14a),  (3.14b)  and  (3.14c)  gives 

Mu.u  55  (3.15a),  ("h)VN/  =  Hv  -t  hj  -2h0  (3.15b), 


F2  Sm2(o(w-trfv) 


F2 


and 


(T)ww  =  i|.a...+-hw — ills 

W  Fz  5m1  (<*Y-  *„) 


(3.15c) 


Equations  (3.15a),  (3.15b)  and  (3.15c)  can  be  expressed  more  conveniently 


as 


(h) 


MM 


+  hw  -  2h 


M 


F2  Si  n  (ofp  — 


(3.16) 


where  M,  N  and  P  have  the  same  meanings  as  were  used  in  equations 


(3.9c)  and  (3.9d). 


The  finite  difference  expression  of  eauation  (3.16)  can  be  sub¬ 
stituted  in  equation  (3.1l)  to  give  an  equation  relating  hQ  to  h^,  h^,  h^,  h”, 
h^  and  h”,  as  follows: 


Z 


+  k. 


( U  >  .  L  11  OL  \  kx  Sin  °<N  •  +  Ky  Co6  oi^j.  Co&ot?  _  ^ 

^nM  +  nM  -  Lno).  --  ■  ,  - - — : — 7 - r - 7 - r  ~  ^  • 

F  Sin  (o(p  Sin(o^Ki  -«W 


M=u3v,w 

Each  term  in  the  summation  may  be  multiplied  by  the  common  factor 


[  F^  sin  (tfp  -of^)  sin 

2  01*  +hM  -2ho) 

M  *u,v,w 


(<*N  -«M)  sin  (o<M  -  c*p)  j  to  qive 

k^  sin  q<n  .SiVio/p  +  co^«xN,  Coze/? 

Sim  (of  h  -  of?') 


-  O.  (3.17) 


Denotina  the  coefficients  C  by 
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^  y,  5»r>  0(\/  5\ n  0(w  +  005  °^v  006  °*  W 

-  — — - - - : - y — : - -  » 

S\r\  C^V-«<W/ 

C  -  kx  Sin  o(wsih  0<u  +•  ky  &>?,  <*„  COS  U a 

Si  n  C6^  w  O 


0  —  5\n  c<u  Sin  o^v  4-  k^  Coao<a  Gbgofy 

Sin  («U-64V) 


equation  (3.17)  can  be  written  as 


Cu.hu+  Cv.h;  +Cw.hi+  Cu..hi+  Cv.hJ  +  CwK-2(0i+Cv+Cw)Ko-O  .  <3-1B) 


Equation  (3.18)  is  the  finite  difference  equation  for  the  head 
at  a  point  inside  a  homogeneous  soil. 


It  is  more  convenient  to  express  the  directions  of  the  Ou,  Ov  and 

Ow  axes  as  the  angle  between  each  and  the  datum  (horizontal)  direction.  If 

these  angles  are  denoted  by  oc^ ,  ot^  and  respectively,  the  coefficients  Cu, 

C  and  C  in  equation  (3.18)  become 
v  w 

Q  _  kx  s'm  (ct y  -0<a)  -(*«*)+  ky  Cosfcy  -ota)  CoS  (o (3.19a) 

Sm  -  o<w) 

Q  _  kx  S>>r>  6»n  (o^u  -otft)-hky  Cos(ul,-O<.o)  Gos(<*i~e(ci)^  (3.19b) 

sin  (o<w  -  <*i) 

£  -  kx  sin  (p<u~o<a)sin(o(v  -o(o)  +  ky  <^(o<u~cOccs(rfv-rfa)  .  (3#19c) 

Sin  w  OL  y) 


u2  axis 


FIGURE  B2.  Boundary  of  two  zones 


Point  on  a  boundary  between  two  zones 

Consider  a  portion  of  a  common  boundary  between  two  zones,  as 
shown  by  the  double  line  in  FIGURE  B2.  The  notation  used  has  the  same  mean¬ 
ing  as  was  used  in  FIGURE  B1  except  that  a  subscript  1  or  2  has  been  added, 
according  to  which  of  the  zones  the  notation  refers.  The  mesh  configuration 
in  each  zone  is  devised  so  that  0,  a  point  on  the  boundary,  is  a  node  common 
to  both  meshes.  One  axis  of  each  mesh  lies  along  the  common  boundary.  In 
FIGURE  B2,  the  Ow  axes  of  both  zones  are  shown  as  the  boundary  axes;  in  general 
any  axis  (Ou,  Ov  or  Ow)  of  one  zone  may  be  matched  with  any  axis  of  the  other 
zone.  Those  axes  not  lying  along  the  boundary  (i.e.,  the  Ou  and  Ov  axes  of 
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FIGURE  ti3 . 


Imaqe  of  zone  1  . 


B8 


both  zones  in  FIGURE  B2 )  do  not,  in  general,  lie  along  the  same  straight 
line  on  either  side  of  the  boundary;  therefore,  although  there  are  six  nodes 
adjacent  to  point  0,  a  finite  difference  expression  similar  to  equation  (3.16) 
cannot  be  stated  directly  for  the  axes  which  are  not  boundary  axes.  Each 
zone  is  considered  separately;  two  ficticious  head  values  at  •image*  points 
(as  shown  for  zone  1  in  FIGURE  B3)  are  introduced  for  each  zone  so  that  fin¬ 
ite  difference  expressions  may  be  stated.  These  ficticious  head  values  are 
treated  as  unknowns  which  are  eliminated  later  by  solving  several  simultaneous 
equations . 


The  boundary  is  shown  in  FIGURE  B3  with  the  image  points  U^  and 
v|  drawn  opposite  Uj*  and  Vj  respectively  by  extending  the  0u^  and  Ov-j  axes 
for  one  mesh  length  in  a  negative  direction.  The  image  points  and  the  fic¬ 
ticious  heads  at  those  points  are  signified  by  the  superscript  *f'. 


By  the  continuity  condition,  the  outward  flow  across  the  curved 
portion  of  the  perimeter  of  the  elemental  semi-circle  on  GH  inside  zone  1  is 
equal  to  the  flow  Q2  into  zone  1  across  the  portion  of  the  boundary  within  GH. 
The  flow  is  expressed  by  integrating  equation  (3.10)  between  (t*wi  -it)  and 


e<w1 ,  that  is, 


w. 

Qi  =  I  <** 

VysM-ir 


'tfw, 

tfWl-TT 


\Vx. 


Cos<*  4  Vyj  Sin  d 


}.  As- 


w, 


nt 


MM,L  2  Jj 


FIGURE  B4 . 


^  Wt 


CoS  Ot  yy  | 


4.  TT  (As)2  y  ^  kjg!  Sl^  (XMi  Sir>  ofp,  +  ^L/l  CoS^M,  CQSo<p, 
2  Sin  (o<ni  ^,n  C°*m  —(*p0 


M3. 20) 


The  velocity  components  at  the  point  J  located  between  G  and  H 
on  the  boundary  (FIGURE  B4 )  are  equal  to  the  velocity  components  at  0  plus 
the  rate  of  change  of  the  velocity  components  along  0w^  multiplied  by  w  (the 
distance  0  J ) ,  i.e., 


velocity  component  at  J  in  -  direction  =  vyi  +  (vy^  )w*wf 
velocity  component  at  J  in  y-|  -  direction  =  Vyi  +  (vyi  )w.w. 


Resolving  these  components  in  a  direction  normal  to  GHf  the  flow  q-j  across 
the  element  JJ*  of  the  length  6w  is  given  by 

%  =  {  (v2:i  +  60w  .  W  )sir>  *w,  -  (vyl  4-  (vyi)w.  w)  coso<wl}iw 


-  (  Vx,  SlVl  <VW,  -  Vyi  cos  0(^1  .  <5w 

+  {  (vjci)w*  <*Wi  ”  (vyi)w*  Co^>  <*  wi}  .  W  .  C^w.  , 


(3.21) 


The  flow  Q2  across  GH  is  obtained  by  integrating  equation  (3.21)  between  -As 


and  +As . 


BIO 


®2L  I  V»l  °^wi  -  Vyi  CO©  ^vg,}.dw 


J.  {(Vxi)w.  Sin  ^  w,  '  (vy«)  w,  Cos  o/w,| .  W  .  dw 


=  ^  { 


V*M  Sin  °<W!  -  V^,  CO SoC^j.&S. 


(3.22) 


Since  =  Q^f  the  right-hand  sides  of  equations  (3.20)  and  (3,22) 


are  equal,  and 


(hW  ^>lf>  .  Sin  (Xpt  -h  k  y  (  QOS  c<N|  Co^t^p) 


,n  C^Ni  SmC°^M«  ”°^p0 


“  O. 


Therefore, 


(k)  kx,Sm  0<M,  SlKl  o(p,  4  ky,  Cos  OtN|CoS  d(p,  _  q 


(3.23) 


M|  =  uovowi 


Sin  <^Mi)*S,r'Co(M,”*o(f*0 


Equation  (3.23)  is  the  continuity  equation  for  flow  inside  the  semi¬ 
circle  on  GH.  This  equation  can  be  written  in  finite  difference  form,  (equation 
3.24),  The  process  of  reasoning  used  for  obtaining  equation  (3.24)  from  equation 
(3.23)  is  similar  to  that  for  obtaining  equation  (3.16)  from  equation  (3.11), 

Qjl.hul‘t"^v»*kYiT'Cvn.hJ|(I-f  Cyt.  hy,  t-CvYthw2.~  2  (Cur+Cy/, +Cwt)ho=:09  (3.24 ) 


where  h'  ,  h,!  r  hr  r  h'  are  the  total  heads  at  U '  ,  V*  W' ,  W' 
u 1  vf  wl  w2  1112 

respectively r 

f  f  •  f  f 

h  .  ,  h  .  ..  are  the  fictional  heads  at  U’  ,  V.  respectively 
ul  vl  11 

and  C  , ,  C  ,r  C  .  are  coefficients  appertaining  to  zone  1  havinq 
ul  vl  wl 

the  same  expansions  as  equations  (3.19a,  b,  c) 


A  similar  equation  exists  for  zone  2;  it  is: 


Bll 


W2.'*"^V2>^y24^'W2‘^w,Z+'QA<2*^U2*'^V2*^^2+^W2.-^W»-2(Ca2+'Cv2+^W^^o®0*  (3.25) 


In  each  zone,  there  are  three  axes;  since  two  axes  are  sufficient 
to  locate  any  point  in  the  plane,  any  distance  along  Ov^  may  be  regarded  as  a 
function  of  ui  and  w-j  only.  Therefore, 

(h)vi  =  (h)u,t(u,)Vl  +  (h)Wl  (w,)v.  •  (3-26) 


But,  by  equations ( 3 .4q  )  and  (3.4h), 


Sin  (dyy/t 

Sin  (o (Wl-olui)  ’ 


(3.27a), 


(3.27b) 


By  expressing  and  in  the  form  of  a  Taylor's  series  expansion,  as  was 

done  for  h^  and  h^  in  equations  (3.12a)  and  (3.12b),  it  may  be  shown  that 

K!u  -  h£  (A.u,)1  aa 

'  2A,u,  6  *  1  'UUU|  . 

Neglecting  terms  containing  AjU^  raised  to  the  power  of  2  or  greater,  finite 
difference  expressions  for  the  first  derivative  of  h  along  Ou-j  f  Qv^  and  Owi  are 


(h)u»  - 

Hu,  -  hi, 

2  A,u, 

5 

(3.28a) 

H 

> 

h  v,  -  h  vi 

2A,v, 

* 

(3.28b) 

1) 

2^ 

^Wl  “  Hw2- 

(3.28c) 

2  AjW, 

Substituting  equations  (3.27a),  (3.27b),  (3.2Ba),  (3.2Bb)  and  (3.28c)  in  equation 


(3.26), 
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Fy,  ~  .  ®'n  6*  Wi  “  Vl)  _  FtWi  ^W2  •  ~  0^ut-tfV(^  .{3.29  ) 

ZA|V,  ZA|U,  Sin^w,-o<ul)  2A|W|  S(o(<Kw,-0(u,) 

By  the  sine  rule  on  triangle  0UjV^t  (FIGURE  B3), 


A,  L» 


I  u  > 


A,v, 


kl  w  i 


s'"'1 *  (c< w, -  V I )  Sin  (<*w,-<*m)  sin  C^VI  -  °<Ui) 

Denoting  each  of  the  above  eaualities  by  F -|  f  and  substituting  for  A|U^, 
and  Aiwi  in  equation  (3.29), 

hyi  -  hyj  -  hu»  ~~  5m(o(\A/i  — c^yi)  _j.  ^vVi  ~  ^ vyq  .  ~Stn{o< ui~ ^vi) 

i  Sin^o^i— t(yi)  (<*vi  Sir\  — ofui) 

Simplifying , 


hu,  -  hj,  -f 


hwi  —  Hj~i  "F  ^v"l  V) 


\A/2 


=r  O. 


(3,30) 


Similarly,  when  considering  zone  2, 


UZ 


-  W\/n  -F  h^-  -F  h  ^ 


U2 


VO,  "  n  VV  \ 


=  O. 


(3.31) 


Any  distance  along  0x^  and  0y^  may  be  regarded  as  a  function  of  u^j 


and  w-j ;  therefore, 


00  XI  —  C^Owi  •  (U  l)  XI  4  (^)wi  *(w,)xi  . 


(3.32) 


But,  by  equations  (3.4c)  and  (3.4d), 


ol  W 


S\n  (o(yiroiUi)  5 


(3.33a),  (W,)^|  -  -  S|n  J  (3-33b> 


Substituting  equations  (3.28a),  (3.28c),  (3.33a)  and  (3.33b)  in 


equation  (3.32), 


■ 
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(h)X)  —  'r>Ul~'^Ul  .  Sin  ^  Wl  ^j_  hyy^  —  h  y/2  .  ~  <*  Ul 

2A,u,  5iir>  (o<W| -tfi*,)  2A<Wi  (<yw,- o(ui) 


Making  the  sine  rule  substitutions, 


hui  -  hfi 


sin  »<wi  ,  hwi  -  hwa  -sin«u.  ,, 

*r  ~  — — : - r  •  — : - — r  y  (3,34) 


^  Xi  2F,  sin(Ww,-o/Vl)*  S^(«wi-oluO  2^1  sin(o<Vl-ofw,)  5»r>^w,-^ui) 

Similarly,  using  equations  (3.4e)  and  (3.4f), 

(h)  -  ^  ut  w  hut  ~Co6  C^WI  hyy/|  —  hyV2_  CoS  <*ui  ^ 

^  2F, sii^ (cVwr o/y,)  uO  2F|  &m  (^yj-c/mi)  Sn(^W|-°(u,) 


The  flow  Q2  across  GH  was  shown  (equation  3.22)  to  be  given  by 

Qi  "  ^  Y^I  5,n  ^WI  -  Coe  oivyi}.  As.  (3.22) 

By  Darcy's  law,  vy1  =  -ky1(h)y1  and  vy1  =  -ky1(h)y1>  therefore 

®2  =  {  ^*1  0*0x1  ^w,  —  ^yi  (fF»)y|  C06  o/W||.As.  (3.36) 

Substituting  equations  (3.34)  and  (3.35)  in  equation  (3.36), 

* 

Qo  =  -2  hufKt,  kxi^^W^l^ytCOS^Wl  _  .  kyiS'in^sin^H-kcf^Sot^Ccfeota,  As  ^  (  g  q7  ) 

2F,  2f,sin(<*vr*ud  5,n  _o(U)^  J 

By  adding  the  expressions  for  Cu1  and  Cy1  (eauations  3.19a  and  3.19b)  it  may  be 


shown  that 


Cm,  +  C 


Ut  ^  W\ 


k-x\  Sir»  ^-yiCpe  d  wi  •  /  >. 

;  ^  ~  ;  ~  ~T  •  SlO  \^Vi—£^U|) 

-o(W|).  S'^  (o^wt  Ut) 


Eauation  (3.37)  may  therefore  be  more  conveniently  written  as 


Q,  = 


=  2  Lj>7h^‘  v(Cu,»Cv,)  +  hwj-hvvi -  Cv,l  As 

(  2r,  5in(cyv, -o^uj)  2F,  Sin  (c^v»-o<ui)  J 

r  {T^ui  “  ^  ifi)- ^Cui4 CVi)-h  (hw,  -hvJ/1).CY, |.(As/A,Vi/,)9 


(3.38) 
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A  further  expression  for  may  be  derived  by  considering  zone  2,  i.e., 

”  ^U2)*  (Cu2  -+  Cv2.  )  +  (hwz  -hwi)Xv2  }  •  (As/A|W2)  9  (3.39) 

Equating  the  riqht-hand  sides  of  equations  (3.38)  and  (3.39),  and  observing 

that  4,-i  -  V2  (-  0W{  =  0W£)  , 


(Cu,^G/,).(hU,-  hi)  +  (Cu2-HCv2)*(hu2*"  ^U2.)  1(3.40) 

*  (Cvi  w  C^va)*  (h  wi  ~  ^W2  )  -  O  .  / 

The  five  equations  (3.24),  (3.25),  (3.30),  (3.31)  and  (3.40)  are 

f  f  f 

mutually  independent,  and  contain  the  four  ficticious  heads  h  h  h 

f 

h  .  These  terms  may  be  eliminated  by  solving  the  five  equations  simultaneously. 


H 


ence 


Cu,  .hy,4  C V,.Kv,+  Cu2^i^+Q2>^^2+'Jfi^^»(hw,+Hw2)-(Cui+Qi4'CwH'^U2+Cy2-hCW2)hoS5  Qj(3. 


41) 


where  Cu1 ,  ,  Cw1 ,  Cu2,  Cv2,  Cw2  are  coefficients  having  the 

expansion  written  in  equations  (3.19a),  (3.19b)  and  (3.19c) 
and  referring  to  the  zone  signified  by  the  subscript. 

Equation  (3.41) is  the  required  finite  difference  equation. 
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APPENDIX  C 


COMPUTER  PROGRAM  :  PROCEDURES  Of  SPECIFIC  PARTS 

Division  of  an  n-  sided  polygon  into  triangles 

An  irreqular  polygon  P^ P£ . . . . Pp  as  shown  in  FIGURE  Cl ,  is 

specified  by  the  coordinates  (x  y  ),  (x  y  ),....,  (x  v  )  listed  in 

1  *1  2  y2  n  n 

clockwise  rotation  round  the  polygon.  Thp  directions  of  each  side  are 
calculated  previously  from  the  coordinates,  paying  due  regard  to  quadrant. 


com- 


A  triangle  is  formed  by  joining  two  alternate  corner  points, 
(i.e.,  P2  to  P^ ,  or  P3  to  P^  ,  etc.).  In  order  that  the  triangle  lie 
pletely  within  the  polygon  the  three  followinn  conditions  must  apply: 

(a)  The  anole  inside  the  polygon  at  the  intermediate  corner  point  is 


less  than  180  degrees.  For  example,  the  triangle  formed  by  joining  P^  to  P^ 

in  FIGURE  Cl  is  inside  the  polygon  since  the  angle  at  P^  is  less  than  180 

degrees;  conversely,  triangle  P  P  P  is  outside  the  polygon  since  the  angle 

5  6  7 

inside  the  polygon  at  P  exceeds  180  deqrees. 

6 

(b)  The  angle  inside  the  polygon  at  the  point  of  the  triangle  which 

bears  the  lowest  subscript  is  greater  than  the  angle  inside  the  triangle  at 

the  same  point.  For  example,  triangle  P  P  P  is  completely  inside  the  polygon 

since  the  polygon  angle  at  P^  (angle  P^P^P^)  is  greater  than  angle  P^P^P^,  but 

triangle  P  P  P  is  not  completely  inside  the  polvqon  since  angle  P  P  P  is  not 
7  8  9  6  7  8 

oreater  than  anqle  P  P  P  . 

7  f  O 

(c)  The  line  joining  alternate  points  is  not  intersected  bv  any  of  the 
sides  of  the  polyqon  which  do  not  pass  throuoh  those  noints.  Trianqle  P^P^P^ 


C2 


C3 
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FIGURE  Cl  Irregular  polygon 


qualifies  on  this  condition,  since  the  line  joining  ^2^4  is  no^  intersected 

by  any  of  the  sides  P^,  P^,  P?PB,  PgP^  P^  ,  but  triangle  P^Pg  is  not 

completely  inside  the  polygon,  since  the  line  P4P6  is  intersected  by  PgP^  and 

P  P 
1  2* 

By  calculating  the  direction  of  the  line  joining  alternate  points 
on  the  polygon,  and  using  the  known  directions  of  the  polygon  sides,  the 
triangle  is  easily  tested  on  conditions  (a)  and  (b).  The  intersection  points 
of  the  line  joining  alternate  points  with  each  of  the  other  polygon  sides, 
produced  if  necessary,  is  calculated;  inspection  of  the  positions  of  these 
points  indicates  whether  condition  (c)  is  satisfied. 

It  may  appear  at  first  that  a  triangle  which  does  not  satisfy 
condition  (b)  could  not  satisfy  condition  (c),  and  therefore  that  condition 
(b)  is  unnecessary.  Although  this  statement  is  true  in  most  cases,  it  is  not 
aenerally  true.  In  FIGURE  C2,  triangle  PpPgP^  which  is  not  completely  with¬ 
in  the  polygon,  fails  to  qualify  on  condition  (b),  yet  satisfies  conditions 


( a  )  and  (c ) . 


C4 


FIGURE  C2 


When  the  date  on  a  polygon  is  ready  to  be  used,  an  inspection 
is  made  to  see  whether  the  polygon  has  only  three  sides.  If  so,  the 
coordinates  of  the  corner  points  are  stored  as  apices  of  a  triangle;  the 
coefficients  Cy ,  C^,  and  Cw  are  calculated,  and  no  further  work  is  reauired 
in  this  part  of  the  program.  If  the  polygon  has  more  than  three  sides,  the 
following  procedure  is  adopted. 

The  points  P2  and  P^  are  joined  and  triangle  P2P3P4  is  inspected 

to  see  if  conditions  (a),  (b)  and  (c)  are  satisfied.  If  not,  P3  and  P5  are 

joined  and  triangle  P^P^P^  is  inspected.  If  this  triangle  does  not  satisfy 

conditions  (a),  (b)  and  (c),  triangle  P^P^P^  i-s  tried,  and  so  on.  The 

geometry  of  any  polygon  dictates  that,  before  triangle  PrP^P^  is  reached,  a 

triangle  satisfying  conditions  (a),  (b)  and  (c)  is  always  found.  The 

coefficients  C  ,  C  and  C  for  this  trianole  are  calculated  and  stored,  to- 
u  v  w 

gether  with  the  coordinates  of  the  triangle  apices.  Calling  this  triangle  by 
the  points  Pr  Pr+-|  pr+2>  where  2  4  r  <  n-1  ,  the  points  Pr+2>  Pr+3»**“*  Pn 
are  re-lettered  pr+1 ,  Pr+2»****»  Pn-1  respectively;  so  that  the  polyaon, 
which  is  reduced  from  n  to  n— 1  sides,  is  lettered 

P^  P^ • • . . pr  and  (re-lettered)  Pr+1  Pr+2 *  *  *  * Pn-1 * 


FIGURE  C3  Part  of  a  flow  region 


The  procedure  is  repeated  until 

reduced  to  three  and  C  ,  C  and 

u  v 

angle . 


the  number  of  sides  of  the  polygon  is 

C  are  calculated  for  the  remainina  tri- 
w 


It  may  be  shown  that, 
polygon  in  FIGURE  Cl  is  divided 


when  the  above  procedure  is  adopted,  the 
into  triangles  in  the  succession  indicated. 


Selection  of  a  triangle  side  which  adjoins  a  given  line 

Before  a  residual  R  can  be  calculated  for  a  node  on  a  side  of  a 
triangle,  the  side  of  the  other  triangle  which  lies  along  the  same  line 
must  be  determined.  It  is  assumed  that  an  inspection  has  already  shown  that 
the  triangle  side  does  not  lie  along  a  constant  potential  face  or  an 


impermeable  boundary  of  the  flow  region 


C6 


Since  the  coordinates  [XU(I),  YU(I)J,  [ XV ( I ) ,  YV(I)],  and 
[XW( I ),  YW(I)]  of  the  apices  of  any  trianqle  I  are  always  listed  in  clock¬ 
wise  rotation,  two  of  the  apices  of  trianqle  J  are  coincident  with  any  two 
of  the  apices  of  an  adjoining  triannle  K  in  one  of  only  three  combinations 


as  follows: 

[XU(J),  YU ( J  )  ] ,  [XV(J),  YV ( J ) J 


or 


[XV(J),  YV(J)],  [XW(J),  YW( J ) ] 
or 

[XW(J),  YW( J ) ],  [ XU ( J  ) ,  YU ( J  )  ] 


may  be 
Dincider 
with , 
respectively 


coincident 

>*  ..  *< 
with , 


[ XW ( K  ) ,  YW(K ) ] ,  fxV(K),  YV(K ) ] 

or 

[XU(K),  YU (K  )  ] ,  [XW(K),  YW(K ) ] 

or 

[ XV (K ) ,  YV(K ) ] ,  [XU(K),  YU (K  )  ] 


If  the  Ow  axis  of  the  given  triangle  lies  along  the  boundary  con¬ 
taining  the  node  under  consideration  and  if  the  number  identifying  the  tri- 
anale  is  5,  the  boundary  passes  through  the  points  [XU(5),  YU ( 5 ) ]  and  [XV<5), 
YV(5)],  (FIGURE  C3)«  By  observinq  the  possible  combinations  for  matching 
apices  of  adjoining  trianales,  the  side  of  the  second  trianqle  which  adjoins 
this  boundary  is  determined  as  described  in  flow  chart  form  in  FIGURE  C4.  For 
the  case  illustrated  in  FIGURE  C3,  the  •'yes”  option  is  taken  at  step  C3,  when 


K  is  6, 


I, 


C7 


FIGURE  C4. 


Flow  chart  for  selectinn  a  trianqlp 
side  which  adjoins  a  given  line. 
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APPENDIX  D 


COMPUTER  PROGRAM  :  FORTRAN  STATEMENTS  AND  TYPICAL  SET  OF  DATA 

FORTRAN  statements 

The  complete  list  of  statements  appears  on  pages  D5  to  D17.  A 
brief  explanation  of  the  more  important  terms  in  the  program  is  given  below. 


A  (  ) 

Array  containing  directions  of  each  side  of  a  polygon  in  radians. 

AA 

Direction  of  minimum  permeability  in  radians. 

AADEG 

Direction  of  minimum  permeability  in  degrees. 

B 

Square  root  of  SM,  divided  by  D. 

BM 

Current  maximum  value  of  B. 

CU(  ) 

CV(  ) 

cw(  ) 

Arrays  containing  coefficients  'C'  for  use 
in  eauations  (3.45)  to  (3.49). 

D 

Maximum  desirable  distance  between  nodes  for  any  zone. 

H(  t  ,  )  Potential  head  at  the  node  signified  by  the  subscripts. 


HA 

Head  acting  on  a  constant  potential  face. 

HAI 

Head  allocated  initially  at  all  nodes  except  those  on  a  constant 
potential  face. 

I 

Triangle  counter. 

IJ 

Counter  used  for  identifying  triangles  when  I  is  in  use. 

IUI  (  )  ] 

I VI  (  ) 
IWI(  )  J 

Arrays  containing  numbers  of  triangles  which  lie  on  impermeable 
[  boundaries. 

IUU(  )  1 

IW(  ) 
IWW(  ) 

Arrays  containing  numbers  of  triangles  which  lie  on  constant 
potential  faces. 

D2 


D3 


11  ( 

12 

14 

15 

17 

J 

JI 

J1 

J2 

J3 

K 

L 

M 

N 

NB 

NIT 

NZ 

P 

PA 

PB 

R 

RF 

RM 

RML 

RT 

R1 

R2 

R3 

5 

5M 


Array  containing  the  lowest  trianqle  number  in  each  zone. 
Total  number  of  triangles  in  the  flow  region. 

^  Counters  used  in  the  format  of  printed  output. 

Iterations  counter. 

Zone  counter. 

Counter  used  for  identifying  triangles  when  I  is  in  use. 
^  Counters. 

^  Counters. 

Number  of  nodes  on  the  side  of  any  triangle. 

Number  of  boundaries  of  any  zone. 

Maximum  number  of  iterations  allowed. 

Number  of  zones  in  the  flow  reqidn. 

TT 

Minimum  permeability  of  any  zone. 

Maximum  permeability  of  any  zone. 

Residual  at  any  node. 

Amount  by  which  the  over-relaxation  factor  exceeds  1. 
Current  maximum  residual  in  current  iteration. 

Maximum  residual  in  the  previous  iteration. 

Tolerable  residual. 

►  Components  of  the  residual. 

Square  of  the  lenoth  of  the  side  of  any  triangle. 

Current  maximum  value  of  5. 


TOL 

TQL 

X  (  ) 
Y(  ) 


^  Tolerances. 

}  Arrays  containing  x-  and  y-  coordinates 
of  corner  points  of  each  zone. 


XC(  ) 
YC(  ) 

XU  (  ) 
YU  (  ) 
XV(  ) 
YV  (  ) 
XW  (  ) 
YW(  ) 


Arrays  containing  x-  and  y-  coordinates 
of  each  node. 


Arrays  containing  x-  and  y-  coordinates 
of  apices  of  each  trianqle. 


XX 

YX 

XI 
Y 1 
X2 
Y2 


^  Coordinates  of  intersection  point  of  two  linps. 

1  Coordinates  locating  a  constant  potential  face 
or  impermeable  boundary. 


SOURCE  STATEMENT 


Sheet  1  of  1 3 


DS 


DIMENSION  I 1( 20 ),X  (  21 ) , Y(  21 )  ,  A(  20)  »XU( 70) , YU (70)  ,XV(70)  ,YV( 70 ) 
DIMENSION  XW(  7  0  ),YW(  70)  ,CU(  70)  ,CV<  70)  ,CW(  70)  ,  XC  (11)  ,YC(  11  ) 
DIMENSION  HI  70,  II,  II  ) 

DIMENSION  IUUI  30) ,  IVVI 30) , IWW( 30) , IUI (30) , I  VI  I  30)  , I Wl  130) 

P-3 • 1415927 
T0L=0. 00001 
TQL=0.0001 

2  HEAD  (5,1000)  N  Z »  R  T »  N  I T 
.  ;  \2  )8, 8, 4 

4  I F I NZ-20 } 9 , 9, 6 
6  WRITE  I  6,  1016) 

8  STOP 

9  1  =  1 

WRITE  (6,1001) 

WRITE  (6,1002) 

BM  =  0 
NGK  =  0 

DO  174  J=  1 , NZ 
Il( J)=I 

SM=0 

READ  (6,1003)  D , P A , P 8 , A ADEG 
I E( P8-PA) 15,  10,  15 
10  WRITE  ( 6,  1004)  J , P  A 
GO  TO  20 

15  WRITE  (6,  1005)  J , P A , P  B , AADEG 
20  AA=AADEG*P/ 180. 

READ  (5,  1006)  X( 1) ,Y(  1 ) 

N6=  2 

25  READ  (5,1006)  X(NB),Y(NB) 

I  F  (  X  (  NB  )  -X  (  D)  35,  30,  35 
30  I  F ( Y( NB )-Y(  1 )) 35,40, 35 

35  NB=NB+ 1 

I  F ( NB-2 1)25,25,  36 

36  WRITE  ( 6, 1017)  J 

37  READ  (5, 1006)  XF, YF 
I  F ( X  F-X (  1)  )37,  38,  37 

38  1F(YF-Y( 1) )37, 39,37 

39  NOK= 1 

GO  TO  174 

40  N B=  NB- 1 

DO  70  K=  1 , NB 
I  F ( X ( K+ 1 ) -X ( K )  ) 50,45,  50 
45  A ( K )  =  P/ 2 . 

I F ( Y ( K  +1  )  -Y ( K )  ) 65,65,  70 
50  A(K)  =  ATAN(  ( Y(K  +  1)-Y(K  )  ) / ( X ( K+ 1 ) - X ( K ) )  ) 

I F ( X ( K  +  1  )  -X ( K )  ) 65, 65,  55 
55  IF(Y(K+1)-Y(K) )60, 70, 70 
60  A ( K )  =  A ( K ) +P 
65  A ( K  )  =  A ( K ) +P 
“70  CONTINUE 
75  L  =  1 

A  B=  A  (  3  ) 

I  F ( NB- 3 ) 148,  148,80 
80  L  =  L  + 1 

I F( L-NB )84, 82, 82  


- 
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82  WRITE  (6*1018)  J 
IM  0  K  =  1 
GU  TO  174 

84  IFt ABSt A(L )-A(  L+l)  )-2.*P  +  T0L)  85,80,80 

85  I F ( A  S  S ( A(L )— At L+i)  )-T0L ) 80, 80,86 

86  IF( AtL)-A(L-H) )89,80,87 

87  I F ( ABS ( A ( L ) ~A{ L + 1 ) -P ) - TOL ) 80 , 80 , 88 

88  I F ( A(L)-A(L+1)-P)91,80,80 

89  IFt ABS( A( L )-A( L+l ) +  P J-TOL )80, 80,90 

90  IF(A(L)-A(L+1)+P)91»80»80 

91  IF(X(L+2)-X(L) )93,92,93 

92  A B=  P/ 2 • 

I F { Y( L+2)-Y(L> ) 96, 96, 97 

93  A 8= AT  AN (  ( Y ( L+2 ) -Y t L ) ) / ( X C L  +  2  ) -X ( L  )  )  ) 

I  F  (  X  (  L  ♦  2  ) -X ( L )  ) 96 , 96,  94 

94  IF(Y(L+2)-Y{L) )95, 97,97 

95  AB= AB+P 

96  AB=AB+P 

9  7  IFt ABSt A ( L  —  1 ) - A ( L )  )-2.*P+T0L)98, 111,111 
9  8  I F ( ABS ( At l-l  )-A(L )  J-TOL ) 1 11 , 1 11 ,99 
99  IFt  At  L-l ) -At  L  >)  102, 11  1, 100 

100  I  Ft  ABS t  At  L-l ) -A (L >-P ) -TOL ) 80,  80,  10  1 

101  I F ( At  L-l ) —  A ( L ) -P ) 104,  80, 111 

102  I F ( ABS ( At  L-l >-A(L )+P  J-TOL ) 80,  80, 103 

103  IFt  At  L-i  )  -  A  (  L  )  +P)  104,  80,  111  _ _ _ 

104  I  Ft  ABS ( At  L-l ) -  A B ) -2 . * P +  TOL )  105, 11 1 ,111 

105  IFt ABSt At L-l )-Ab)-TOL  )  111,111, 106 

106  I  Ft  At  L-l )-AB) 109, 11 1, 107 

107  IFt  ABSt  At  L-l J-AB-P )-TOL ) 80, 80, 108 

108  IFt At L-l )-AB-P )  111, 80, 80 

109  _ I F(ABS(A(L-1)-AB+P J-TOL ) 80, 80 , 1 10 _ 

110  IFt  At  L- 1 J-AB+P )  111,80,80 

111  M=  L  +  3 

IFtL-NB+l ) 113,  112,  112 

112  M  =  2 

GO  TO  115 

113  IF(M-NB) 116, 116,114 

1 14  M=  l 

115  IFt M-L+2) 116, 1 16, 148 

116  XY=(X{MJ-XtM+l))*(Y(L+2)-Y(L))+(X(L+2)-X(L))*(Y(M+l)-Y(M)) 
IFtXY ) 1 17, 146,  1  17 

117  X(M)*YtM+l)-X(M+l)*Y(M) 

DL=  X(L  +  2)*Y(L ) — X ( L ) *  Y ( L  +  2 ) 

Y  X  = ( DM* ( Y { L+2)-Y< L) )+DL*(Y(M+l)-Y(M) ) )/XY 
XX=(DM*tX(L+2)-X(L) )+DL*(X(M+l)-X(M) ) )/XY 
IFtABStYtMJ— YtM+l) )-ABS(X(M)-X(M+l)J )123,118,118 

118  IFt ABSt YX-YtM)  J-TQL ) 128,  128, 1  19 

119  IF(ABStYX-YtM*l))-TGL)l28,128,120 

120  IF(YX-YtM) ) 1 2 1  ,  128,  122  _  _ _ 

12 1  I  Ft  YX-Y (M+l )  ) 146, 128,  128 

122  IF(YX-Y(M+1) )128,  128,  146 

123  I  Ft  ABS ( XX-X ( M)  ) -TQL )  1 2 8 ,  1 28 , 1 24 

124  IFtABStXX-XtM+l  )  )-TQL  )128, 128,125 

125  IF(XX-XtM)  )  126,  128,  127 

126  IF(XX-X  (M-t-1  )  J146,  128,  128  


' 
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127  I F(XX-X (M+l  )  ) 128,  128,  146 

128  iHABS(Y{  L  )  -Y  (  L +2  )  )-AQS(X(L)-X  ( L+2 ) )  ) 134 » 129  ,  1 2  9 

129  I F ( ABS( YX-Y(L) ) — T QL ) 80, 80, 130 

130  IF(A6S(YX-Y(L+2))-TQL >80,80,131 

131  I F  <  YX -Y ( L  )  )  132, 80,  133 

132  IF(YX-Y(L+2j >14  6,  80,  80 

133  I E( Y  -Y ( L+2  )  )80 ,00,  146 

134  1  Ft  »BS ( XX-X ( L )  ) -TQL ) 80, 80,  135 

135  I F { AOS { X  X-X { L  +  2  )  ) - TQL ) 80, 80, 136 

136  I F { XX-X ( L  )  ) 137,  80,  1  38 

137  IF(XX-X(L+2) 1146,80,80 

138  I F ( XX-X ( L+4  )  )£0 , 80,  146 
146  M-M+l 

IF(M-L)  115,113,113 

148  S  = ( X ( L  +  2  )-X l L )  )  **2  +  { Y { L  +2 ) - Y ( L  )  )  *  *  2 
I F( SM-S  )  150,  152 , 152 
150  S M=  S 

152  S=(X(L+1)-X(L>  >**2  +  (  Y  {  L  +  1  )-Y <  L )  ) *  *  2 
IF( SM-S )  154,  156,156 
154  S M=  S 

156  S= ( X ( 1  +  2 )-X ( 1+ 1  )  ) **2  +(Y(L+2)-Y(L+l)  )  ** 2 
IF( SM-S  )  158,  160,160 
158  SM=S 

160  XU (  I  ) =X { L  +  1  ) 

YU( I )=Y(L+1 ) 

XV  i I ) =X i L  +2 ) 

YV(  I ) - Y (  L  +  2  ) 

'  X  W  (  I  )  =  X  (  L  ) 

YW (  I  )  =Y  CL) 

S  N=  P A*S  I N  (  A  (  L  )  -  AA  )  *  S  I N  (  A  {  L  +  1 )  -  A A  ) 
CS=P8»C0S(A(L)-AA)*C0S(A(l+1)-AA) 

CU( I )  =  ( SN+CS ) / S  IN ( A ( L +  ) 

SN=PA*S  IN(  A(  L  +  1  )-AA  )  *S  IN  (  AB-AA  ) 

CS=PB*COS  i  At  L  +  1  )-AA  )  *COS(  AB-AA  ) 

CV( I )  =  ( SN+CS >/S  INI A(L  +  1)-A8) 

SN=PA#S IN { AB-AA  )*S IN ( A( L )-AA ) 

CS=P8'-C0S  (  AB-AA  )*COS(  A(  L  )-AA  ) _ 

CW (  I  1  =  1 SN+CS ) / S  IN (  AB— A ( L ) ) 

1=1  +  1 

I F ;  1-70  )  161, 161,167 

161  i r ( NB- 3 ) 170, 170,162 

162  A(L)=AB 

164  I  F  (  L-NB  +  2  )  166,  166,  168 _ _ _ 

166  X(L+1 )  =  X ( L  +2  ) 

Y ( L  +  l  )  = Y ( L  +2  ) 

A ( L  +  1 )  =  A ( L  +2 ) 
l  =  L  +  l 

GO  TO  164 

168  X(NB)=X(NB+1)  _ 

Y ( NB ) -Y ( NB+  1  ) 

N6=  NB- 1 
GO  TO  75 

167  I F ( NB-3 )  169,  169, 171 

169  I F  C J-NZ ) 171,  170, 170 
171  WRITE  (  6,  1019)  
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STOP 

170  12=1-1 

B=SQR  T ( SM )/D 
I F ( BM-8 )  172,  174, 174 
172  8M=  8 
174  CONTINUE 

I F ( NGK ) 8, 175,8 
1  75  N  -  B M  +  2  . 

IF(N-li ) 178, 178,176 
176  N= 1  1 

1 78  DO  184  1=1,12 

DO  182  K  1=  1 ,  N 
N I  =  N-K 1  +  1 
DO  180  K2=  1 , N I 
180  H( I ,K1,K2 )=-999. 

182  CONTINUE 

184  CONTINUE 
WHITE  (6,1007) 

WHITE  (6,1008) 

I  3=0 

J  1=  1 
J  2=  1 
03=1 
14  =  0 
H  A  A  =  0 

186  HEAD  (5,1009)  X  l , V 1 , X 2 , Y 2 , HA 

I F ( ABS ( HA +999. ) -TQL ) 1 89 , 1 89 , 1 87 

187  14=14+1 

H  A  A  =  H  A  A  +  H  A 
I F ( X  2-X  1  )  196,  1  88,  196 

188  iF( Y2-V  i  )  L96,  190,  19o 
169  WRITE  ( 6, 1020) 

STOP 

1 90  J 1 1= J 1—1 
J21  =  J2-  1 
J  3  1  =  J  3-  1 
TOT= 14- 1 
HAI=HAA/TOT 
Du  195  1=1,12 
DO  193  K  1=  1 , N 
N I  =  N— K 1  +  1 

DO  191  K2=  1 »  N I 

I F ( ABS ( H(  I , Kl, K2) +9  99.  ) - TQL  J 1 85 , 1 8 5 , 1 9 1 

18 5  H (  I , K 1 , K2 )  =  HA  I 

191  CONTINUE 
193  CONTINUE 
195  CONTINUE 

WRITE  (6,1011) 

WRITE  (6, 1008) _ 

13=1 
J  1=1 
J  2=  1 
J  3=  1 

192  READ  (5,1012)  X1,Y1,X2,Y2 
I  F ( X 2-X 1 ) 196,  194,  196  


SOURCE  STATEMENT 


SHpr-t  OT  1  .1 


194  I F ( Y  2-Y  1  }  196,246,  196 
196  WRITE  (6,1010)  X1,Y1,X2,Y2 
.  197  DO  198  1=1,  12 

IF(X1-XU(  I  )  )  20  6 , 200 , 206 

200  I  F { Y i-YU ( I )  )206, 202,206 

202  I  F ( X2-X V (  I  >  )  19  8,204,  198 

204  IF(Y2-YV(  I))  198,222,  198 

2  06  : ~ ( X 1  — X  V (  I  )  ) 21 4, 208, 2 14“ 

208  1 F ( Y  1-Y V (  I  )  )214,210,214 

210  I  F ( X2-XW (  I  )  )  198,212, 198 

212  I F ( Y 2 - Y W (  I  )  )  198,230,  198 
2  14  I  F ( X 1 —X  W (  I)  )  198,216, 198 
216  I  F ( Y  1- Y W (  i)  )  198,218, 198 
218  IF(X2-XU(  I  ))  198,220,  198 
220  I  F ( Y2-YU(  I  )  )  198 ,238, 198 

198  CONTINUE 

199  WRITE  (6,1020) 

STOP 

222  I  F ( Jl-30)223,223,  171 

223  I F(  I  3  )228, 224,  228 

224  DO  226  K1=1,N 
NI=N-K1+1 

226  *H  (  I  ,  K  1  ,  N  I  )  =  HA 
IWWt J1)=I 
J1=J1+1 
GO  TO  243 
228  .1  W  I  (  J1)  =  I 
J1=J1+1 
GO  TO  192 

230  IF(J2-30)231»231, 171 

231  I F (  13)236,  232,  236 

232  00  234  K1=1,N 
234  H(  I  ,K1,  1  )=HA 

IUU ( J2 ) = I 
J2= J2+ 1 
GO  TO  243 
236  I U  I  (  J  2  )  =  I 
J2=J2+L 
GO  TO  192 

238  I  F ( J3— 30  )239,239,  171 

239  I  F (  13)244,240,  244 

240  DO  242  K2=1,N 

242  H ( I , 1 ,K2)=HA 
I  V  V  (  J  3  )  =  I 

J  3= J  3  + 1 

243  DO  233  1=1,12 

I F ( X 1— X  U (  i  )  )  20 3, 201, 203 

201  I  F ( Y  1-YU (  I )  )  20 3, 207, 203 

203  I  F ( X2-XU (  I  )  )209,205,209 

205  I  F ( Y 2-YU (  I )  )  209 , 207 , 209 
207  H {  I  , 1,N  )  =  HA 

GO  TO  233 

209  I  F ( X  1  —  X V (  I  )  )213,211,213 

211  I  F ( Y 1 - Y  V (  I)  )213,217,213 

213  I  F ( X 2-X V (  I )  >219,215,219 
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215  IF(Y2-YV(  I )  )219f2i7,219  _ 

2  L  7  H  (  I  ,  N  ,  1  )=HA 
GO  TO  233 

2  19  IHX1-XW1  I)  )  22  5, 2  2  1,2  25 
221  I  F  i  Yi-YW(  I  )  1225,229, 225 
225  IF(X2-XW(  I  )  1233,227,233 
227  I  F ( Y2-YVH  I  )  1233,229,233 
229  H(  I  #  ft  1  1  =  H  A 
233  C  TINUE 

GO  i U  186  t 
244  I  V  I (  J  3  1  =  I 
J  3=  J  3+  1 
GO  TO  192 
246  j 12= J  1- I 
J22  =  J2- 1 
J  32  =  J  3-  1 
I  7  =  0 
R  F  =  0  •  8 

_ 247  RM=0 

DU  586  1=1,12 
DO  584  K  1  = 1 , N 
N  I  =  N— K 1  +  1 
-DO  582  K 2  = 1 »  N I 
I  3  =  0 

248  I  F ( K1+K2-N-1  1250, 326,  250 
250  I  F ( K  L -  1  1252,290,252 
262  .  IF (K2-1 1254,256,254 

£54  R  1  =  CU  (  I  1  *  (  H  (  I,K1+1,K2)+H(  I  ,  K  1- 1 ,  K2  )  1 
R2  =  CV (  I  )  * ( H(  I ,  K  1,  K2+  1  1  +h(  1  ,  K  1  ,  K2-  1  1  1 
R3=CW (  I  )*( ri(  I, K  1+1, K 2- 1 ) +H(  1  ,K  1-  1  ,K2+  1  )  ) 

R  =(R1+R2+R3)/(2.*(CU( I 1  +  C  V  ( I  1 +C  W ( I 1  ) 1  -  H  (  I ,K1  ,K2) 

GO  TO  578 
256  J2=0 
258  J 2=  J 2  +  1 

IF( J21-J2 1262,  260,260 
260  I F(  IUU( J2  1- I  1258,  574,  258 

_ 262  J2  =  0 _  _ 

264  J 2= J 2  + 1 

IF( J22-J21270,  266,266 
266  I F(  IU  I  ( J2  )- I  1264,  268,  264 
268  R  1  =  CU  (  I  1  *  (  H  (  I,K1—  1,  1  1  +  H  (  I,K1+1,1)  1/2. 

R2=CV (  I  )*H(  I, K 1,2)+CW(  I)*H(I,K1-1,2) 

GO  TO  366 _ _  _ _ 

270  NK1=N~K1  . 

DO  272  I J=l,  12 
I  F ( XV ( I  ) -XU (  IJ  )  1276,273,276 
273  I F ( YV ( I  )-YU (  I J  1  1276, 2 74,  2  76 
2  74  I  F( XW (  I  )  —  X W {  IJ  1  127  2,2  75,2  72 

275  I  F { YW {  I  l-YW (  I J  1  1272, 234,  272  _ _ _ 

2  7  6  T  F  (  Ml  1-  X  V  (  I J  )  1280, 27  7, 2  oO 

277  I F { YV (  I  l-YV (  I J  1  1280, 2 78,  280 

278  IF(XW(  I  1  —  XU (  IJ  1  1272,279, 272 

279  I F( YW (  I  1  — YU (  IJ  1  1272, 286, 272 

280  I  F( XV (  I  )  —  X W {  IJ  1  1272, 281, 272 

281  I  F  (  YV  (  I  )-YW(  IJ  1  1272,282,272 
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232  I F ( X W (  I  )  — X V (  IJ  )  )272»283»272 
28  3  I F ( Y W (  I  )  — Y V (  I J  )  )272, 238, 272 
272  CONTINUE 
GU  TO  199 

2  84  R1=(CU(  I  )+CV(  IJ  ))*(H(  I,K1-1,1)+H(  IJ*l,Kl+l))/2. 

K2  =  CV {  I  ) *H(  I ,K 1 , 2 ) +CW(  I J ) *  H  C I J , 2f K 1-1 ) 

R3=CW (  I  )*H(  I,K1— 1,2)+CU(  IJ)*H(IJ,2*K1) 

GU  TO  364 

286  R 1= ( CU  (  I  )  +  CW( J J  ) ) *  (  HI  I,Kl-l,l)+H(IJ,Kl+l,NKl))/2. 

R2  =  CV  (  I  )  *H(  I  ,  (s  1 ,2  )  +CU(  IJ)*H(IJ,K1-1,NK1  +  1) 

R3=CW( I ) *H( I,K 1-1 i 2 )+CV( IJ )*H( I J ,K1,NKI ) 

GU  TO  364 

288  R 1=  C  CU (  I )+Cg(  IJ  )  )*( HI  I ,K1-1, 1 )+H( I J,NK1 , 1) ) /2. 
R2=CV(  I  )*H(  I , K 1 , 2 ) +C V l I J )*H(  I J,NK1  +  1 ,2) 

R  3  =  C  W (  I  ) *  H {  IjK 1-1, 2  )  +  C  W  ( IJ )  *  H  ( IJ,NK1,2) 

GO  10  364 

290  I F ( K2-1 )292, 368,292 
292  J  3  =  0 
294  J  3=  J 3+ 1 

IFIJ31-J 3)298, 296,296 
296  I F (  IVV( J3)-I )294,  574,  294 
298  J  3=  0 
300  - J3= J3+1 

I F ( J  32- J  3 ) 306,  302,  302 
302  I F (  IV  I ( J3 )- I )300,  304, 300 
304  R 1  =  C V (  I  )*{ H(  I,  1 ,K2+1)+H(  I ,  1 , K 2- 1 ) ) /2 . 

R2=CU (  I  )*HC  1 , 2*  K2 ) +CW (  I  )*H{  I , 2,K2-  1  ) 

GU  TU  366 
306  NK2  =  N-K 2 

DO  308  I J  =  1 ,  I ? 

I F ( X  W (  I  ) —  X  V (  IJ  )  )  3 1 2  *  309,312 

309  IF(YW( I )-YV( IJ ) )312,310,312 

310  I F< XU(  I  )-XU(  IJ  )  )308, 311, 308 

311  IF< YU< I  >  — YU (  IJ  )  )308, 320, 308 

312  I  F ( X w (  I  )-XW(  IJ  )  )  3  16, 313, 316 

3  j  3  I F ( YW (  I  )  -Y  W (  i  J  )  )  3 1 6 , 3  1 4 , 3 1 6 
314  I F ( XU (  I  ) -X  V (  {J  )  )  30  8,  315,  308 
3^5  I F ( YU (  I  )  —  Y  V (  I J  )  ) 308, 322, 308 
3  16  I F  <  X W (  I  ) -X  U (  IJ  )  ) 308, 3 1 7,  308 

I  F  (  Y  w!  (  I  )  -  Y  U  (  IJ  )  )  3  0  8  ,  318,  308 

318  I  F ( XU (  I  )-X W (  I J  )  )  30  8 ,  3  19,  30  8 

319  I  F ( YU (  I  )  - Y  W {  IJ  )  )  308, 324, 308 
308  CONTINUE 

GO  TO  199 

320  R 1 = ( C V (  I  ) +CW (  I J  ) )*( H(  I , l,K2+l )+H( I J , NK2  +  2 , K2- 1 )  ) /2. 
R2=CW(  I  >*H{  1,2, K2-1 )+CU( I J ) *H ( I J , NK2 , K2 ) 

R3=CU(  I  )*H{  I,2,K2)+CV(  IJ  )*H( I J,NK2  +  1,K2-1) 

GO  TO  364 

3  22  R 1= ( CV (  I  )  +CU( IJ  )  ) * { H ( I , 1 , K2+1 ) +H ( I J , K2- 1 , 1 ) )/2. 

~  R2=CW ( rT*HT  iT27K2-l)+CV( IJ)*H{ IJ,K2,2) 

R3=CU ( I  )*H( I,2,K2)+CW( IJ)*H(IJ,K2-1,2) 

GO  TO  364 

324  R  1= ( CV {  I  )  +  CV (  IJ  )  ) * ( H (  l,l,K2+l)+H(IJ,l,NK2+2)  )/2. 

\  R2=CW (  I  )  *H(  1 , 2, K2-1 )+CW(IJ)*H(IJ,2,NK2) 

R3=CU (  I  )  *H(  I , 2 , K2 ) +CU (  I J  )  *H(  I  J , 2 , N K 2  +  1  ) 
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GO  TO  364 

326  I  F ( K 1- 1 )  328 ,436,328 
328  IFCK2-1  )330,  504,330 
330  J  1  =  0 
332  J  i  =  J  1  + 1 

I  F ( J 1 1- J  1  )  336,  334,  334 
334  I  F  (  fW  (  )i)-I  )  33  2 ,  574  ,  332 
3  36  J  I  ~ 

338  Jl=Jl+i 

IFU12-JL)  34^  340,  340 
340  I F{  rw I { J 1 )- I ) 338, 342, 338 

342  R  1  =  CW  (  I  )  *(  H(  I,  K  1-1 ,  K2  +  1  )  +  H(  I  ,  K  1  + 1  ,  K2- 1  )')  /2. 

R2=CU( I )*H( I,K 1-1,K2)+CV( I )*H( I ,K1 ,K2-1 ) 

|  GO  TO  366 
344  DU  346  IJ=  1,  12 

I F{ XU (  I  )-XW (  IJ  )  )350, 34  7, 350 
347  IF(VU( I  )-YW(  lJ  )  )350, 348,350 
J48  I F ( XV (  I  )-XV (  IJ  )  )346,  349, 346 

349  I F ( Y  V (  I  J  —  Y  V (  IJ  )  )  346 ,  3  5  8,  346 

350  IF(XU( I  )-XU(  IJ  )  )354,  351,  354 

351  I F ( YU (  I  )-YU(  IJ  )  >354,352,354 

352  I  Ft  XV (  I  )-XW(  IJ  )  )346, 353, 346 
3  53  •  I F( YV (  I)-YW(  IJ  )  )346,  3  6  0,  346 
354  I F( XU(  I  )-XV(  IJ  )  )346, 355, 346 
3  55  I F ( YU (  I  )  —  Y  V (  I  J  )  )346,3  56,34  6 
3  56  I F{ XV (  I  )-XU(  IJ  )  )  346 ,  3  5  7,346 
35  7  I F  <  YV (  I  )  — YU (  I J  )  )346,  3  62, 346 
346  CONTINUE 

GU  TO  199 

358  R  1= ( C W (  I  > *CU(  I J  )  ) * { H(  I»Kl+l»K2-l ) +  H ( I J , K 1 - 1 ,1 )  )/2. 

R2= CU( I ) *H( I, < 1-1, K2)+CV(IJ)»H(IJ,K1, 2) 

R  3  =  C  V  t  I  )  *H<  I,K.1,K2-1)+CW(  I  J  )  »  H  (  I  J  ,  K  1-1 , 2  ) 

GU  TO  364 

360  R  1  =  (  C  W  (  I  )+CV(  IJ))*(H(I,Kl+l,K2-l)+H(IJ,l,K2+l))/2. 

R2=CU (  I  )  *H(  I  ,< 1-1 , K2 ) +CW(  I J )*H(  I J, 2 , K2-  1  ) 

R3=CV(  I  )  *  H (  I,<1»K2-1)+CU(  I  J  )  *  ri  (  I J , 2  »  K2  5 
GO  TO  364 

362  Rl  = ( CW (  I) +  CW( I J  )  )* ( H(  I , K 1+ 1 , K2- 1 J +H ( I J , K2  +  1 , Kl-1 )  )/2. 
R2=CU(I)*Fi{ I,K1-1,K2)+CU(IJ)*H( IJ,K2-1,K1) 

R3=CV(  I  )*H(  I ,K 1 ,K2-1 )+CV(  IJ )*H{  IJ,K2,K1-1) 

364  R  ={Rl+R2  +  R3)/lCU(  I  )+CV(  I  ) +C  w  { I  )+CU(IJ)+CV(lJ)+.CW(IJ)  ) 

R  =  R  -H(I,K1,K2) 

GO  TO  578 

366  R  =(R1+R2)/(CUI  I  )+CV(  I  )+CW( I  )  )-H( I ,K1,K2) 

GO  TO  578 

368  I F ( 13)402,370,402 
370  R 1=0 
R2=0 

I  J=  I _ _  _ _ 

372  J  I  =  I J 
J  2  =  0 

374  J2= J2+ 1 

I F ( J21-J2J378, 376,376 
376  I F (  I UU ( J  2 ) — J I ) 374,  574,  374 
3  78  R 1  =  R 1 +CU ( J I )  *  H  ( JI,2, 1  )  +  C  V  {  J  I )*H( J I , 1 , 2  ) 
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R2=R2+CU ( J I ) +  C  V  ( J  I  ) 

J  2  =  0 

330  J 2  =  J 2+  l 

I F ( J  2  2- J  2 ) 384,  382, 382 
382  I Ft IU 1 ( J2 )-J I) 380, 572, 380 
384  DO  386  I  J=  1 ,  12 

I E ( XV { J  I  )  -XU (  IJ  )  )  39  1,  387,391 
387  IF ( YV( J  I  > -YU (  I J  )  ) 391, 388, 391 
3  88  I F ( X W ( J  I  )  — X  W (  IJ  ) ) 386,  389, 386 

389  IF(YW( J i )-YW( ! J  )  )386,  390,386 

390  I F (  I J - I  )  3 72 ,  576,372 

391  I F( XV ( J I )-XV( I J  )  >396,  392,396 

392  I  Ft  YVt  J  I  )-YV( I J  )  ) 396,  39  3, 396 

393  I  F ( X  W ( J  I  )  — X  U (  l  J  )  ) 3  8  6 ,  394,  38  6 
39  4  I r ( Y W ( J  I  )-YU(  IJ  )) 386,  395 , 386 

395  I F { I J- I ) 440 , 576,440 

396  I  F ( XV { J  I  )-XW(  I J  )  ) 386,  397, 386 

397  IF(YV( J I )-YW( I J  )  ) 386, 398, 386 

398  I F ( X  W  <  J  I  )-XV(  IJ  )  ) 3  8  6 ,  399,  38  6 
3  99  I  Ft YW( J  I  )  — Y V ( I J  )  ) 386, 400, 386 
400  I F (  I  J- I  >508,576,508 

386  CONTINUE 
.GO  TO  199 
402  IJ=I 
404  J  1=  I  J 
J  3=  0 

406  J  3  =  J  3+  I 

IF(J31-J3)4I0, 408,408 
408  I F {  I V V ( J3 )-J I ) 406, 574, 406 
410  I F (  13 >414, 412, 414 

412  R i  =  R 1 +  CV ( J I )  *H  (  J I  , 1 , 2  )  +CU ( J I )*H( JI ,2, 1 ) 
R2  =  R2  +  CV( J  I )  +CU  <  J  I  ) 

GO  TO  416 
414  I  3  =  0 
416  J  3  =  0 
418  J  3=  J  3+ 1 

IF( J32-J3 >422, 420, 420  _ _ 

"420  I F ( IVI ( J  3  )-J I ) 4  18, 576,418 
422  DO  424  I J  =  l,  12 

1  F  (  X  W  (  J  I  )  -XV (  I J  )  1428,425,42  8 

425  I F { YW ( J  I  )  — Y  V  <  IJ  )  1428,426,428 

426  IFtXUt JI  )-XU( I  J  )  ) 424, 427, 424 

427  IF(YU(J I  )-YU( IJ  )  ) 4  2  4 ,  540,424 _ _ 

428  1  F ( XW ( J I )-XW( IJ  )  >432,429,432 

429  IFtYWtJ  I  )  — Y W (  l J  )  >432,430,432 

430  I F ( XU ( J  I  )-XV (  I J  > > 424,  431 , 424 

431  I F { YU( J I  )-YV(  I J  > >424,404,424 

432  IFtXWtJ  I  >-XU(  I J  )  >424,433,424 

433  I F ( YW ( J I ) —YU (  I J  )  >424,434,424  _ _ 

434  I  Ft  XU ( Jl ) —  X W { I J  )  >  424,435,424 

435  IF(YU(J  I  J-YvJ  (  IJ  )  >424, 472,424 
424  CONTINUE 

GO  TO  199 

436  I F ( 13)470,438, 470 
438  R  1  =  0  
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R2  =  0 
I  J*  I 

440  J I = I J 
J  3  =  0 

442  J3=J3+l 

I  F ( J  3  1  - J  3 ) 446»  4  44, 444 
444  I F (  I V V ( J3)-J I) 442, 574,442 

446  R  1  =  R  1  +CV ( J I  }*H( Jl,  1,N-1)+CW(JI  )*H( JI  , 2,N-1) 
K2=R2+CV ( J I  )  +CW ( J I  ) 

J  3  =  0  i 

448  J  3  =  J  3+  1 

I  F ( J32-J3 )452,  450, 450 
450  I  F(  IV  I  ( J3  )-J I ) 448,  572, 448 
452  DO  454  I J=  1 ,  12 

I F ( X  W { J I  ) -XV (  I J  ) ) 459, 4  5  5, 459 

455  I F ( YWt J I )-YV( I J ) )459, 456,459 

456  I F ( XU ( J I )-XU( I J  )  >454,457,454 

457  If (YU(JI  ) — YU  C  IJ  ) )454,  458,454 

458  I F{ I J-I )440, 576,440 

4  59  I F ( XW { J I )-XW ( I J  )  ) 464, 460, 464 

460  IF(YW(J I ) — YW( IJ ) ) 464 , 461,464 

461  If ( XU ( J I )-XV( I J ) )454, 462,454 
462 • If (YU( J I ) -Y V ( I J ) >454, 463,454 

463  If (IJ-I )508, 576,508 

464  If ( X W ( J I ) -X  U { I J  ) 1454,465,454 

465  If (YW( J  I  ) -YU (  I J  )  ) 454, 466 , 454 

466  If ( XU( J  I  )-XW(  I J  )  )454,  467,454 

467  I  f ( YU ( J I  )-YW( I J  )  )454, 468,454 

468  If (  IJ-I  )  372, 57  6, 372 
454  CONTINUE 

GO  TO  L 99 
470  .  I J= I 
472  J  1=  I J 
J  1  =  0. 

474  Ji=JI+I 

IF.(Jil-Jl)478,  476,476 
476  IF(  IWW( J  I  )-J I ) 474, 574, 474 
478  If (  13)482,480,  482 

480  R  i  =  R I +CW ( J I  )*H ( J  I  , 2*N- 1 ) +CV l J I ) *H{  JI  , 1 ,N-1 ) 
R2=R2+CW( J I ) +CV ( J I ) 

GO  TO  484 
482  13=0 
484  J 1=0 
486  Jl= Ji+i 

If  (  JI2-JD490,  488,488 
488  If (  IW  I(J  i  )-J I ) 486,  576,486 
490  DO  492  I  J=  1 ,  12 

If (XU( JI )-XW(  I J  )  )  496, 49  3,496 

493  If ( YU( J  I  )-YW( I J  )  )496,  494,496  _ 

49  4  If(XVTTlT=XV(  I  J  )  )  492 , 49  5,492 
49  5  I F ( YV { J  I  )  -Y V (  I J  )  ) 49  2, 404,492 

496  I  F ( XU ( J  I  ) -XU ( I J  )  )  500, 49  7, 500 

497  I f ( YU ( J  I  )-YU ( I J  )  ) 500, 498 , 500 

498  I  F ( XV ( J  I  )-XW (  I J  )  ) 492, 499, 492 
4  99  If ( YV( J  I  )-YW(  I J  )  ) 49  2, 4  72, 492 


' 
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500  I F ( XU ( J  I  )-XV(  IJ  )) 492, 501,492 

501  i  F ( YU { J I )  -  YV  ( IJ  )  )492* 502,492 
50  2  I  P  (  X  V  (  J I  )-XU( IJ  )  )  4  9  2  *  50  3,492 
50  3  IF(YV< J  I  )-YU(  IJ  )  )49  2,  540,4  92 
492  CONTINUE 

GO  TO  199 

504  I F (  1  3)538, 506,  533 
506  R 1=0 

R2  =  0 

r  i  —  r 

m  vJ 

508  JI=IJ 
J  1  =  0 

510  Ji=Jl+i 

IF(  Jll-Ji ) 5 14,  512,512 
5.12  IE  (  .  WW  (  J  1  )-J  t )  5  10,  5  74, 510 

5  14  RL=R1+CW( JI  )*H( Jl,N-l,2)+CU( JI ) * H ( JI  , N- 1  ,  1 ) 
R2=R2+CW(  J  I  )  t'CU  ( J  I  ) 

J  1  =  0  <r 

516  J  1  =  J 1  + 1 

I  F ( J12-J  1  ) 520,  5  18,  518 
518  I  F  (  IW  I  {  J  1  )-J  i!)  5  16,  5  72, 516 
520  DO  522  IJ=1,  I-fc 

-IF  (XU  (  J  I  > -X  W  ( I  J  )  )  52  7,  52  3, 52  7 

523  IF( YU( J I )-YW(  i  J  )) 52  7, 524,  52  7 

524  I F ( X V  <  J  I)-XV(  I J  )  )522, 525,522 
5  25  I  F ( Y V ( J  I  )  -Y  V ( I J  ))522, 526,  522 

526  I  F { IJ-I )508, 576,508 

527  I  F ( XU ( J  I  )-XU( I J  )  )  532, 528, 532 
5  28  I F ( YU ( J  I  )  -YU ( I J  ) ) 5 32 , 529 , 5 32 

529  I  F ( XV ( J I )-XW( I J  ) ) 522, 530, 522 

530  I F ( YV( J I )-YW( I J  )  )522, 531, 522  J _ 

531  I  F (  IJ-I  )372, 576,372 

532  I  F ( XU ( J  I  ) —XV ( I J  ))  522, 533, 522 
3  33  I  F ( YU ( J  I  )— YV  C I J  )  )  522,  534, 522 
334  I  F ( X  V ( J  I  ) -XU ( I J  ))522, 53  5,522 

535  IF( YV( J  I  ) —YU (  I J  )  )522,  536,522 

536  I F (  IJ-I  )440, 576,440  _ 

p~22  CONTINUE 

GO  TO  199 
538  I J  = I 
540  J I = I J 
J2  =  0 

542  J2= J2+ 1 

'  I  F  ( J21-J2 )  546*  544,  544 
544  I  F (  IUU( J2  )-J I ) 542, 574, 542 
546  I F ( 13)550, 548, 550 

548  R 1  =  R 1 +CU ( J  I )*H(JI*N-1,1)+CW(JI)*H(JI,N-1,2) 
R2  =  R2  +  CU ( J I ) +C  W ( J I  ) 

GO  TO  552 
550  13=0 
552  J 2  =  0 
554  J2= J2+ 1 

I  F ( J  2 2- J2)558,  5  56,556 
556  I F( IU I ( J2 )-J I ) 554, 576, 554 
558  DO  560  I J  =  1 ,  12  
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I F( XV ( J I ) - X U {  I J  )} 564, 561 , 564 

561  IF(YV(JI)-YU(IJ ))564, 562,564 

562  I F ( XW { J I)-XW( I J  )  ) 560, 563, 560 

563  1 F ( YWi J I  )-YW{  I J  ))560,472,560 

564  I F ( XV  C  J I ) —XV  C I J  ))568,565,568 

565  1F( YV( J I )-YV( I J  )  ) 568,  566, 568 

566  I  F { XW { J  I  )-XU ( I J  )  )  560, 567,560 

567  I  F ( YW( J i  ) — Y  U  t  [J  )  )  560, 540,5 60 
5  68  I  F { X V ( J  I  ) — X W {  I J  ))  560,  5  65, 560 
5  69  I r ( Y  V ( J I  ) - Y  W {  I J  )) 560,  570,  560 

570  I F ( X  W ( J 1 >  —  XV {  IJ  } )560, 571,560 

571  I F ( Y W ( J I  )-YV(  I  J  )  ) 560 , 404 ,  5 60 

5  60  CONTINUE  _ _ _  _ 

GU  TO  199 

572  13=1 

GO  TO  248 
574  R=0 

GO  TO  578 

576  R  =  R  1/  R  2-H (  I  ,K1,K2) 

578  I  F( RM-ABS( R>  )580,  582, 582 
580  RM=ABS(R) 

582  H(  I,Ki,K2)=H( I,K1,K2)+R*( l.Q+RF ) 
584  -CONI INU6 
586  CONTINUE 

I F ( 17-1 )592, 592,588 
588  I F ( RM-RML ) 592, 592, 590 
590  RF=RF  *0 . 8 

592  RML=RM 
17=17+1 

I F { 17-NIT ) 593, 594, 594 

593  I F( RT-RM )247,594, 594  _ 

594  WR I f £  (6,  102 2)  17, RM 

J=1 

14=2 

15=50 

DO  604  1=1, 12 

I  F ( J-NZ  )  595, 595 , 598  _ _ 

595  IF (  I  1 ( J  >-I  ) 598,  596,  598 
.596  I  F (  15-49 ) 599, 599,  597 

597  WRITE  (6, 1013)  J, 14 
14=14+1 

WRITE  (6,1014; 

15=5 

'  GO "10  601  . 

599  WR  I  TE  (6,  1021  ) 

WRITE  (6,1024)  J 
WRITE  (6,1014) 

I  5=15  +  5 

601  JC= J  _ 

J  =  J  +  1 

598  DO  602  K1=1,N 
N I  =  N-K 1  +  i 

DO  600  K2=  1 »  N I 
E  K  1  =K  1-1 
E  K2  =  K  2-  1  
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EN=N- 1  _  J__ _ _  __  _ _ 

XC(K2)=XW( I )+( XV( I )— X W ( I ) )*EK1/EN+ ( XU ( I )-XW( I ))*EK2/EN 
600  Y  C  (  K  2  )  =  Y  W  (  I  )  +  (  Y  V  (  I  )-YW(  I  )  )*EK‘1/EN  +  (  YU  {  I  )-YW(  I ) ) *EK2/EN 
WRITE  (6,1015)  (XC(K2) ,YC(K2)  ,H( I ,Ki,K2)  ,K2=1,NI ) 

I F ( N I —8  )  606  7  60 6  »  608 
606  I F( NI-4 )612» 612,610 
bQ8  15=15  +  1 
610  15=15+1 
612  15=15+1 

IF(  15-56)602  fb 02,614 
614  I  F* (  N — K  1  )  6 1 6 »  6 1 1> »  62  0 

616  I F (  I  2- I  )  604 ,  604,617 

617  IF( J-NZ ) 6 18, 618 ,620 

618  I F (  I  1 ( J  ) - I - 1 ) 6  2  0 , 604,620 
620  WRITE  ( 6,  102  3] f JC,  14 

14=14+1  ^ 

WRITE  (6,1014) 

15=5 

602  CONTINUE 

I  F (  I  5-5)622,604,622 
622  WRITE  (6,  102  1 ; 

I  5=  I  5+  L 
604  CONTINUE 
GO  TO  2 

1 QOO  FORMAT! 7X13, FI 0.3, I  10) 

100 10  FORMAT (  1 H  14  1X4  H  ZONE  8X  7HM  IN  IMUM8X7HMAXI  MUV-5X9HDI  RECTI  G.\|4X 
18HPAGE  1) 

1002  FORMAT (  1H  <+3X46HN0  PERMEABILITY  PERMEABILITY  OF  MIN  PERM/) 

1003  FORMAT (  IXF9.2, 2E10. 2, F  10.2 ) 

1004  FORMAT  (  1H  42'Xl  3, 1PE15.3,  3X9HI  SO  TROPIC  ) 

1005  FORMAT (  IH  42X I  3, 1P2E 15. 3, OPF 14.2 ) 

1006  FORMAT (  1XF9.2, F  10. 2  ) 

1007  FORMAT ( 1H545X40HL0CAT IONS  OF  FACES  OF  CONSTANT  POTENTIAL/) 

1008  FORMAT (  LH  47X 1 HX9X ,  1HY  14X ,  1HX9X , 1H Y/ ) 

1009  FORMAT ( IXF9.2,4F10.2) 

1010  FORMAT (  1H  40X2 F  10 . 2 , 5X2F  10 . 2 ) 

1011  FORMAT (  1H748X35FLUCAT IONS  OF  IMPERMEABLE  bOUNDAR I E  S / ) 

1012  FORMAT (  1XF9.2,  3F10.2 ) 

1013  FORMAT (  IH162X4HZGNE I  3, 4 2X4 H PAGE  14/ ) 

10i4OFGRMAT(lH  9X  LHX 9X IHY 7X4HHEAU 1 IX 1HX9X 1HY7X4HHE AD  1 1 Xl HX 9X 1H Y7 X4HHE AO 
11 IX  1HX9X  LHY7X4HHEA0/  ) 

1015  FORMAT ( IH  F12.2,2F10.2,F13.2,2F10.2,F13.2,?F10.2,F13.2,2F10.2) 

1016  FORMAT ( 1H158X14HTOO  MANY  ZONES) 

1017  FURMATi  iH  50X28HT00  MANY  BOUNDARIES  FOR  ZONE  I  3) 

1018  FORMAT  (  IH  54X2  0  PiDAT  A  INCORRECT,  Z0NEI3) 

1019  FORMAT (  1 H  54X2 3HSECT I  ON  TOO  COMPLICATED) 

1020  FORMAT (  IH  58X1 4HDAT A  INCORRECT) 

1021  FORMAT (  IH  ) 

1022  FORMAT ( IH« 41X2 2HL ARGEST  RESIDUAL  AFTERI4,15H  ITERATIONS  WASF8.3) _ 

1023  FORMAT (  IH 156X4HZUNE I  3 ,  12H  (  CONTINUED ) 36X4HPAGEI 4/) 

1024  FORMAT (  IH  62X4 HZON E I  3 / ) 

END 
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Typical  set  pf  data 

The  first  card  of  each  set  contains  the  number  of  zones,  the 
tolerable  residual,  and  the  allowable  number  of  iterations.  The  data  for 
each  zone  follows  next,  and  appears  in  batches  of  cards,  each  batch  con¬ 
taining  all  the  data  for  one  zone.  The  first  card  of  each  batch  contains 
the  maximum  desirable  distance  between  nodes,  the  minimum  and  maximum 
permeabilities  (in  exponential  form),  and  the  direction  of  minimum  permea¬ 
bility  measured  counter-clockwise  from  the  horizontal.  The  remaining  cards 
each  contain  the  coordinates  of  the  corner  points  of  each  zone  in  clockwise 
order;  starting  and  finishing  with  the  same  point. 

Following  the  zone  data,  one  card  appears  for  every  constant  poten¬ 
tial  face.  Each  card  contains  the  coordinates  of  the  points  between  which  the 
face  lies,  i.e.  XI,  Y1 ,  X2,  Y2,  followed  by  the  head  HA  acting  on  the  face. 

A  dummy  card  consisting  of  zeros  follows  the  constant  potential  face  data. 

The  data  relating  to  the  impermeable  boundaries  follows  next,  each 
card  containing  the  coordinates  of  the  points  between  which  the  boundary  lies. 
A  dummy  card  consisting  of  zeros  completes  the  set  of  data. 

A  further  dummy  card  containing  a  zero  follows  the  final  set  of 
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FIGURE  D1 .  Typical  set  of  data 


A  typical  set  of  data  is  shown  in  FIGURE  D1  .  This  data  relates 
to  a  flow  region  in  which  there  are  3  zones.  The  minimum  and  maximum 
permeabilities  of  two  zones  are  6  x  10”^  and  24  x  10“^  cm/sec,  and  0  x  10“^ 
and  32  x  10“^  cm/sec  respectively  and  they  are  horizontally  stratified,  i.e., 
the  direction  of  the  minimum  permeability  is  90  degrees.  The  third  zone  is 
isotropic,  and  has  a  permeability  value  of  64  x  10“^  cm/sec.  The  tolerable 
residual  is  0.01  and  the  maximum  number  of  iterations  permitted  is  120. 
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TYPICAL  SET  OF  OUTPUT 

As  each  set  of  data  is  analysed,  the  output  is  printed  in 
tabular  form.  A  typical  example  of  output  is  presented  on  paces  E3 
to  E7.  The  information  therein  is  plotted  in  FIGURE  5.3fc0in  the  text. 

The  first  paae  of  each  set  of  output  consists  of  the  data  des- 
cribinn  the  flow  reqion,  the  number  of  iterations  nerformed,  and  the 
larnest  residual  in  the  final,  iteration. 

Succeedinn  oaqes  contain  tabulations1  of  the  potential  head  values 
at  every  node  in  the  flow  reaion.  There  are  four  tables  across  each  paqe, 
and  each  table  consists  of  three  columns.  The  numbers  in  the  first  and 
second  columns  are  the  x-  and  y-  coordinates  of  a  node;  the  potential 
head  at  the  node  is  qiven  in  the  third  column. 
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ZONE 

MINIMUM 

MAXI  MUM 

DIRECTION 

PAGE 

NO 

PERMEABIL I TY 

PERMEABILITY 

OF  MIN  PERM 

Of  OUTPUT 

1 

6 . 000E-06 

2 • 400E-05 

90.00 

2 

8. 000E-06 

3. 200E-05 

90.00 

3 

6.400E-05 

I SOTROPIC 

LOCATIONS  OF  FACES  OF  CONSTANT  POTENTIAL 


X  Y 


X  Y 


0.  100,00 

150.00  100.00 


50.00  100.00 

200.00  100.00 


LOCATIONS  OF  IMPERMEABLE  BOUNDARIES 


X 

Y 

X 

Y 

• 

200.00 

100.00 

200.00 

75.00 

200.00 

75.00 

200. 00 

50.00 

200.00 

50.00 

200.00 

0. 

200.00 

0. 

0. 

0. 

0. 

0. 

0. 

50.00 

0. 

50.00 

0. 

75.00 

0. 

75.00 

0. 

100.00 

50.00 

100.00 

150.00 

100.00 

LARGEST  RESIDUAL  AFTER  120  ITERATIONS  WAS 


0.010 
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THE  FINITE  DIFFERENCE  EQUATION 
AT  AN  APEX  OF  SEVERAL  TRIANGLES 


F  1 


APPENDIX  F 


DERIVATION  OF  THE  FINITE  DIFFERENCE  EQUATION 
FOR  THE  POINT  AT  AN  APEX  OF  SEVERAL  TRIANGLES 

A  derivation  of  equation  (3.43)  in  the  text  was  obtained  after 
this  thesis  was  written.  This  derivation  is  presented  below. 

Consider  a  polygon  ABCDE  enclosing  the  apex  0  of  several  triangles, 
as  shown  in  FIGURE  FI ,  where  the  Ou  and  Ov  axes  of  each  triangle  form  the 
boundaries  and  the  Ow  axes  are  not  shown.  The  points  Uj,  U^ ,  U^,  U^,  U^  are 
nodes  adjacent  to  the  point  0  along  the  Ou  axes  of  triangles  1,  2,  3,  4,  5  res¬ 
pectively.  A,  B,  C,  D,  E  are  situated  at  the  mid-points  of  OUj,  OU^,  OU^,  OU^, 
OU^  respectively  so  that  the  linear  dimensions  of  the  polygon  ABCDE  are  one- 
half  of  those  of  the  polygon 

As  in  SECTION  3.4,  use  is  made  of  the  concept  that  any  distance  along 
the  principal  permeability  directions  Ox^  and  Oy t  and  also  the  direction  Ov^  , 
in  trianqle  1  can  be  regarded  as  a  function  of  distances  along  the  Ou^  and  Ov^ 
axes  to  give 


OOxt  ~ 

(^)  Ul  • 

(“0 

•  xt 

+ 

Mwi  • 

(FI 

a ) 

(H)y|  - 

G3)U|  . 

6*,; 

lyi 

■+* 

(^)wi  * 

(FI 

b ) 

and 

(h)vi  = 

Mui  • 

>VI 

•h 

63)  W1  • 

(w,)V;  • 

(FI 

c ) 

Substituting 

expressions 

for  (u.| 

^  x  1  » 

U 

'1  >x1 

,  (W^yT,  (Ui)v1, 

(w1  >v1 

obtained  from 

equations 

(3.4c  to 

h) 

in 

equations  (Fla,  b,  c), 

F2 


. 


TRIANGLE  2 


v,  ax i6 


TCiANGLE  3 


vs  axis 


TRIANGLE:  -4. 


TRIANGLE  f 


F  3 


TRIANGLE  5 


V4  axis 


FIGURE  FI,  The  apex  of  several  triangles. 


(h)*,!*  0"Wf: 


n  o( 


w  I 


S'n  C^wi-^ui) 


0) 


Wl 


Sin  c* 


Ul 


6m  (d  vvi  — °fu\} 


Hi “”^,7  + 

Ai)y,  =  (h)iW.  a)n  - 

Sin  (o(  w«  ”**u0 


^  .  CoS  cX  u  i 

v '  V  Wl  *  ;  /  - 1_3 — r  ? 

Sin  0*  wi  -  o<  wi_) 


w,  .  s'n  0*w -■**,)  . 

Sin  (of*,  -of*,) 


The  term  (h)w^  is  eliminated  from  equations  (F2a,  b,  c),  giving 


(h)x,-(H) 


sm  oC  vi 


(h)y,=“Wui 


CoS  cXvi 
sin  (oty,  -«ut) 


-  (H). 

+  (h)vi 


Sin  c><ui 


U’  5irt  fa?,  -o^u,)  V^V|  s,n  (oCv,-o«u.i)  ’ 


Cos 


Oi. 


Sin  (oiVi  -  o4t J 


(F2a) 

(F2b) 

(F2c) 


(F3a) 

(F3b ) 


Assuming  that  the  velocity  components  vx1  and  vy-|  in  the  principal 
permeability  directions  are  uniform  along  the  length  AB  in  FIGURE  FI,  the 
flow  across  AB  is  found  by  resolving  the  components  vx1  and  vyi 


in  a  direction 


• 

F4 


perpendicular  to  AB  and  multiplying  the  result  by  the  length  AB  (=  ^.Aiw-j); 
that  is , 

flow  across  AB  =  ( V^.,  «*Wl  -  Vy(  cos  *  Wl)  .  ,J2.  A  |W,  , 

which,  by  Darcy's  law,  is  equal  to 

—  ^in  <*WI  ”  (h)yi 

where,  in  this  expression,  ( h )  x  ^  and  (h)y-j  are  the  partial  deriva¬ 
tives  of  h  at  points  along  AB. 


Substituting  the  expansions  for  (h)y^  and  ( h )  ^ ,( equations  F3a,  b), 
in  the  above  expression,  the  flow  across  AB  is  equal  to 


S\o  ofy,  61  n  q<wi 

(o<v, 


o(  qi  SlO  Q^vV ) 

5in  (oiy,  -o<ui) 


^yi  (h\ii  <LOS>  ^vi  Cc>s  °*wt  —  *lu  c&S  (A|Wl))t(F4 ) 

5tr\  Cct-\j\  —  <^ui)  £o<V(  ~o<Ui)  J  Z 

where,  in  this  expression,  (h)y^  and  (h)v-j  are  the  partial  deriva¬ 
tives  of  h  at  points  along  AB.  These  terms  are  given  approximately  by  the 
finite  difference  expressions 

(k)ul  *  -hui  ~ho  ,  (F5a),  (h)v|  =  hi  -ho  (F5b) 

A,n,  A,v, 

By  the  sine  rule  on  triangle  OUjlV, 

At  u,  ^  Aiv,  -  A,w, _ . 

Sin  (p<W|-^v,)  Sin  (o<w,-^uO  Si" 

Using  these  relationships  and  the  fact  that  h^  =  h^t  equations  (F5a,b)  are 
altered  to 

(h)u,  =  ^Ut  ~  ^ o  .  (o^yt  -  txQ 

A,W,  (^wi-c^vi) 


* 


(F6a) 


. 


F5 


and 


(h)v,  -  ^ut.  ~  .  3*^  (o(yi 

A|W,  (o<vVi  *"  oi.  ui  ) 


(F6b) 


Substituting  equations  (F6a,b)  in  equation  (F4),  and  rearranging, 
Flow  across  AB  = 


{• 


—  H ~ (  C(y(  S>H  o(yV)  t  k^|  CerS  C(y^  CQ5 Oi^j\ 

A,w,  sin  (o^yy,  —oiv,") 

+  &l,a  <*u|S»nc<v\/l  +■  k-ui  CoSo(W(  I  *  (A,w<) 

I  2 


A,wt  5«to  (e<Wi  ~  <*ui) 

Cui  ^ lm  +■  ^vi  ^ui  ""  (Cu,  ■+’  Cvj)  In, 


(F7) 


where  ,  C  ^  are  coefficients  appertaining  to  zone  1  having  the 
same  expansions  as  equations  (3.19a,  b). 

Continuity  requires  that  the  net  flow  across  sides  of  the  polygon 
ABCDE  be  zero;  hence,  by  summing  expressions  similar  to  equation  (F7), 


Cu,4-Cyg  ^|+  Cu2  +  Cv,  hg^  Gus+Cvah^  Cut+Cvs  h  C<J5+Cv4  hn 
Z  Z  Z  tZ  2 


US 


-{ 


Cut  +  Cu2'+’  Cug-t*Cu4.  'f’CuS  +  Cyt  -+CW+-Ctf5  +  CvH+Cysl  K  =o 

2  J 


This  equation  can  be  stated  in  the  form  of  equation  (3.43)  for  cases  in  which 
there  are  n  triangles  with  apices  at  0. 


(Cm+Cvr>)  hg,  *+*  (Cui+^V»)^U2  + - - +  (C  un  +  ^vrvi)  kn  1(3. 

“  (c ui  +  Cv«  -♦■Cuz  +  Cvt.  +  - 4-  Cun*^  ^vn)^o  = 


43) 


